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Abstract 


Fibre Bragg grating sensors have emerged as a simple, inexpensive, accurate, sensitive 
and reliable platform, a viable alternative to the traditional bulkier optical sensor 
platforms. In this work we present an extensive theoretical analysis of the tilted fibre 
Bragg grating sensor (TFBG) with a particular focus on its polarization-dependent 
properties. 

We have developed a highly efficient computer model capable of providing the 
full characterization of the TFBG device in less then 3 minutes for a given state of 
incident light polarization. As a result, the polarization-dependent spectral response, 
the held distribution at the sensor surface as well as the fine structure of particular 
resonances have become accessible for theoretical analysis. As a part of this computer 
model we have developed a blazingly fast full-vector complex mode solver, capable of 
handling cylindrical waveguides of an arbitrary complex refractive index profile. 

Along with the theoretical study we have investigated optical properties of the 
TFBG sensor with application to polarisation-resolved sensing. We proposed a new 
method of the TFBG data analysis based on tracking the grating transmission spectra 
along its principle axes, which were extracted from the Jones matrix. 

In this work we also propose a new method of enhancing the TFBG sensor refrac- 
tometric sensitivity limits, based on resonant coupling between the TFBG structure 
resonances and the local resonances of nanoparticles deposited on the sensor surface. 
The 3.5-fold increase in the TFBG sensor sensitivity was observed experimentally. 
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Chapter 1 


Introduction 


In the presented work we study the tilted fibre Bragg grating (TFBG) sensor and 
present new methods of improving its sensitivity. The TFBG sensor has proven to 
be a simple, inexpensive, accurate, sensitive, and reliable platform [6], with a variety 
of applications, including chemical and biological sensing [7, 8]. 

The TFBG sensor is based on a standard telecommunication fibre with a tilted 
grating inscribed inside the core. The telecommunication fibre is an optical cylindrical 
waveguide consisting of a high refractive index core surrounded by a cladding layer 
with a lower refractive index. The light propagating in the fibre core is confined by 
the total internal reflection phenomenon. 

In our research the 1-cm-long TFBGs were inscribed in the hydrogen-loaded core 
of a standard telecom single mode optical fibre (Corning SMF-28) using the phase 
mask technique and intense ultraviolet light (pulsed KrF excimer laser at 193 nm or 
248 nm). The TFBG sensor coated with small spherical particles is shown schemati¬ 
cally in Figure 1.1. 

The inscribed grating planes are slightly tilted relative to the fibre cross-section, 
which allows to couple the forward-propagating light from the fibre’s core to the 
backward-propagating cladding modes [9, 10]. More than a thousand propagating 
cladding modes are usually excited. Each cladding mode can be viewed as a propa¬ 
gating electromagnetic wave, with a unique wave vector and field distribution, thus 
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Perturbation in the refractive index 



Cladding 


Figure 1.1: Schematic representation of the TFBG sensor coated with sensitivity 
enhancing layer of nanoparticles. 

interacting differently with the outside medium. 

The coupling between the core and the cladding modes, mediated by the grating, 
produces discrete narrow attenuation bands in the transmitted spectrum (the “reso¬ 
nances”) and since the cladding modes are guided by the fibre-surrounding medium 
interface, these resonances are extremely sensitive to the refractive index of the sur¬ 
rounding medium [11], The energy coupled to a particular cladding mode as well as 
the corresponding propagation constants can be precisely measured by acquiring the 
transmission spectrum of the grating. The spectrum is shown in Figure 1.2, consisting 
of a set of narrow resonances sensitive to the refractive index change of environment. 

The TFBG sensor has a unique set of advantages, including the possibility to 
isolate temperature effects and develop compact and robust refractometers with a 
wide operating range [11, 12]. In addition, the simultaneous probing of the external 
medium at various wavelengths with modes having unique polarization properties and 
incidence angles is possible. 

The sensitivity of the TFBG sensor can be further improved if modes of the 
sensor can be coupled to an external resonant system, such as a nano-scale coating, 
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Figure 1.2: A typical spectrum of a 10 degree TFBG sensor. 

with properties sensitive to the surrounding medium. The cladding modes excited 
by the TFBG structure have a non-zero evanescent field at the cladding boundary, 
therefore the light can tunnel outside the fibre into a coating layer. If the phase 
matching condition is met the resonant coupling between the cladding modes and 
resonances of the coating layer can occur. When some modes (among the large number 
of cladding modes accessible) are phase matched to the surface plasmon resonance 
(SPR) of the outer surface of the metal coating, the light energy incident from these 
cladding modes is efficiently coupled to the surface plasmon-polariton oscillations [7]. 
The energy coupled to the SPR causes a drastic change in the TFBG transmission 
spectrum at wavelengths corresponding to the resonant modes. The principle of 
operation of the TFBG SPR sensor is analogous to the well-known Kretschmann- 
Raether setup [13, 14], except that the cladding modes use wavelength and angle 
scanning to probe for the SPR phase matching condition (because each cladding 
mode couples at a different wavelength and “strikes” the cladding boundary at a 
different incident angle). 

Finally, a further increase of the selectivity of the modes coupling to SPR is 
provided by controlling the light polarization. This discovery led the way to a strong 
increase in the accuracy with which we can follow SPR shifts associated with small 
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refractive index changes of the outer medium surrounding the fibre by using the 
polarization-dependent loss (PDL) of the TFBG transmission [15]. In particular, 
it was found that narrowband resonances (100 pm spectral bandwidth (BW)) with 
refractometric sensitivities (S) of 350 nm/RIU (refractive index units) showed up 
in the PDL spectrum. These resonances have some of the highest figures of merit 
(S/BW = 3500 RIU -1 ) for SPR sensors of any kind [16]. 

Several resonances in the transmission spectrum are shown in Figure 1.3. The 
resonances were continuously observed during the process of nano-scale metal coating 
deposition. The complicated structure of the resonances, with strong dependence on 
the him thickness, is proportional to the time of the deposition is clearly seen. 



t(sec) t(sec) 

Figure 1.3: Evolution of the TFBG spectral response during the silver nanoparticle 
deposition followed by the continues him formation. 

Recently it was discovered that the sensitivity of TFBG sensors can be improved 
by coating its surface with randomly oriented silver nanowires. Such a coating would 
allow for a large number of cladding modes to interact with the deposited nanowires, 
exciting localized surface plasmon resonances (LSPR) when the appropriate phase 
matching condition is met. The sensitivity of the TFBG resonances to the external 
refractive index change was increased by a factor of 3.5 even though the surface 
coverage of the nanowires was less than 14%. 
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In the presented work we investigate methods of improving the TFBG sensor 
sensitivity by coating its surface with various types of nano-scale films, including 
metal films and nanoparticle coatings. We also discuss the polarization properties of 
the sensor and present the developed highly sensitive measurement technique. 

However, the main objective of our work was to simulate the behaviour of the 
TFBG sensor. The sensor is an interesting object with more than a thousand inter¬ 
acting modes and a complex polarization-dependent response, as can be seen from 
ures 1.2 and 1.3. The simplicity of the spectral response hides many non-trivial 
physical effects that we might wish to decode. Our goal here was to understand the 
TFBG sensor behaviour, interpret the polarization-dependent spectral response and 
create a computer model which would allow us to predict the outcome of experiments. 
Moreover, the process of improving the sensor sensitivity required the rather tedious 
experimental work of scanning the parameter space by probing particles with different 
geometrical shapes, such as spheres, cubes, cages and wires, varying the nanoparti¬ 
cle’s size and the deposition densities. A theoretical guidance which would allow 
the prediction of an optimal set of parameters, such as the him material, thickness, 
composition and morphology is highly desirable. 



The thesis Organization 


This work is organized in the following manner. Chapter 1 provides an introduc¬ 
tion to the problem and motivation. 

In Chapter 2 we explain how we developed a full-vector complex mode solver 
which we applied to circularly symmetric optical waveguides in order to compute 
scalar and vectorial modes. The exact and approximate solutions were compared and 
the orthogonality of the modes was verified. An interesting analogy between TE and 
TM modes in slab waveguides and modes in cylindrical structures was established. 
The analysis of a weakly-guided approximation and the exact solution allowed us to 
build a connection between the scalar and vectorial equations. The split in eigenvalues 
was discussed and analyzed thoroughly. To our knowledge, this particular analysis 
has not been reported previously. 

Although a number of commercial mode solver software packages is available, the 
software interfaces usually impose certain limitations on the process automation. For 
example, to simulate a TFBG spectral response more than a thousand interacting 
modes should be considered. These modes should be computed, often an various 
frequencies, stored and be available for further processing. If the optimal parameters 
are searched, a number of separate computations must be done sequentially. In 
addition, in order to study a particular property of the sensor only the modes with 
this particular property are required, thus these modes can be targeted discriminately 
without computation of the remaining modes. Thus, by targeting only specific modes 
the simulation speed can be increased dramatically. 

Chapter 3 is based upon solutions obtained in Chapter 2. The modes obtained 
for the non perturbed case were used as basis functions to solve the problem of tilted 
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Bragg grating structure. The polarization-dependent effects were investigated. The 
results of the simulations were compared with the experimental measurements. The 
emphasis was set on polarization dependency of the coupling coefficients, which were 
computed for all possible angles of incident core mode polarization, to our knowledge 
this result has not been previously reported. 

Chapter 4 provides a complete experimental characterization of the TFBG struc¬ 
ture. The experimental measurement techniques with application to polarization- 
based sensing were developed. 

In Chapter 5 we start our discussion on the possible enhancement of TFBG sensor 
sensitivity by coating its surface with nano-scale films. We review optical properties 
of various materials, metals in particular. 

Chapter 6 provides a detailed overview of nanoparticle optical properties, as well 
as several simulation methods, in particular the Mie theory. 

In Chapter 7 we present our ideas on enhancing the sensitivity of TFBG sensors, 
and conduct a search for parameters that would provide the optimal nanoparticle- 
based coating. 

Chapter 8 presents experimental results from various nano-scale film coatings and 
their associated sensitivity enhancements. 

Lastly, Chapter 9 is a summary of the work findings. 




Chapter 2 


A full-vector complex mode solver 
for circularly symmetric optical 



2.1 Introduction 

In this chapter we present a simple yet efficient and exact approach for fast and ac¬ 
curate modeling of waveguides with cylindrical or elliptical symmetry, with arbitrary 
real and imaginary refractive index profiles. 

The analytical solution for the simplest case of a cylindrical dielectric waveguide, 
as well as modes classification and field plots were presented by Snitzer in 1961 [17]. 
Shortly after, direct numerical integration techniques were presented, in a noticeable 
paper published by Dil and Blok [18] on a numerical solution of four coupled Maxwell’s 
equation in radial parabolic dielectric waveguides. 

It is well known that in the case of waveguide with a small refractive index contrast 
the weakly guided approximation can be used, allowing the system of Maxwell’s to be 
reduced to the Helmholtz equation, which can more easily be solved [19, 20]. Here we 
consider the general problem of full vectorial mode solutions of waveguides with high 
refractive index contrasts and a non-zero imaginary part of the refractive index, i.e. 
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we consider an arbitrary complex permittivity profiles. 

Several methods have been proposed earlier. Usually the waveguide profile is 
subdivided into a number of piecewise homogeneous concentric layers, and the system 
of Maxwell’s equations are solved in each layer. The solutions are next connected with 
the help of boundary conditions for the tangential electric and magnetic fields at the 
layer boundaries. 

In cylindrical coordinates the solutions for each layer can be represented in terms 
of Bessel and modified Bessel functions, and then connected through a 4x4 transfer 
matrix, which incorporates the boundary conditions. This is the so-called transfer 
matrix method proposed by Yeh and Lindgren [21]. The modes propagation constant 
are obtained by finding the roots of a polynomial, over Bessel functions, with a degree 
proportional to the number of layers. This method becomes increasingly prohibitive 
with regards to computational resources as well as to numerical stability with the 
increase in number of layers. The transfer matrix method becomes even more com¬ 
plicated if the permittivity posses an imaginary part or the modes are leaky. In such 
a case the roots have to be searched on the complex plane. 

Although its limitations, the transfer matrix method can be successfully applied 
to waveguide structures consisting only of a few uniform layers [22, 23, 24] or to 
structures consisting of infinitely many alternating uniform layers, the so-called Bragg 
fibres [25]. The computational complexity of the matrix method can be reduced 
by replacing Bessel functions with their asymptotic expressions, or if the periodic 
cylindrical layers are approximated as planar Bragg stacks [26]. 

Another possibility is to use the pseudospectral method, also called the spec¬ 
tral collocation method, where a solution is searched in terms of sine functions [27], 
Chebyshev-Lagrange [28, 29], Laguerre-Gauss [30, 31], Her mite-Gauss [32] interpo¬ 
lating polynomials, or in terms of some other suitable basis functions. The pseu¬ 
dospectral method can be implemented either for the entire domain or separate basis 
functions can be chosen for each uniform layer. In each case the basis functions are 
chosen with regards to boundary conditions. If the problem is solved for the entire 
domain the approximation of mode fields might require hundreds or even thousands of 
functions, and the convergence is generally a problem [33]. The multidomain method 
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is in a way similar to the transfer matrix method. The dielectric interface conditions 
should be fulfilled, thus the functions between different layers are again connected 
through the transfer matrix at the boundaries [29]. In both cases the system of alge¬ 
braic equations results in numerical eigenvalue problem, with eigenvectors consisting 
of the expansion coefficients for a particular mode. A simpler but more resource 
hungry method was proposed in [34, 35] where instead of the global functions for 
each homogeneous region a finite different method was applied, and then the trans¬ 
verse electric or magnetic fields were matched at the radial index discontinuities. The 
comparison between the transfer matrix method and the pseudospectral method was 
conducted in [33, 29] and definitely favours the pseudospectral method. 

Yet another possibility is to use a standard finite difference (FDM) or finite element 
(FEM) method for the entire 2D waveguide profile, with the Yee’s cell [36] adapted for 
the cylindrical symmetry [37, 38, 39]. A number of commercial simulation software 
based on the finite difference and finite element methods are available. This approach 
does not require understand of engineering or physical aspects of the problem. A given 
problem is viewed as a “blackbox” accepting the input data e.g. geometrical shape, an 
operational wavelength, etc. and returning the required numerical results. However, 
the use of commercial software imposes certain limitations, for example it might be 
challenging to automate the process of obtaining a series of solutions for a varying 
parameter. It should also be noted that the FDM and FEM method are the slowest 
and the most resource hungry, especially if used for 2D or 3D problems. 

Here we present a simple yet efficient and fast numerical method. First the system 
of Maxwell’s equations was reduced to only 2 coupled ordinary differential equations 
for the electric held. The variation of the dielectric permittivity was incorporated in 
the equation such that equations become suitable for numerical integration through 
the entire domain, in contrary to the piecewise homogeneous methods described 
above. Next the equations were transformed into a system of algebraic equations 
with the help of a finite difference method. Finally the problem was reduced to find¬ 
ing eigenvalues and eigenvectors of a five-diagonal matrix. The eigenvalue problem 
was effectively solved with the standard iterative method commonly applied to large 
sparse linear systems. The proposed method requires less grid points if the refractive 
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index discontinuity is approximated with a smooth function at an interval of about 
1 nm at the position of discontinuity. 
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2.2 The solutions for a cylindrical waveguide 

We start with Maxwell’s equations written in the MKS system of units, for a 
source free region [40]: 

V x H{r, t ) = e -^p-d t E(r, t ), (2.1a) 

V x E(r, t ) = d t H(r ., t), (2.1b) 

c 


V-H(r,t)= 0, (2.2a) 

V ■ (e(r)E(r, t)^j = 0. (2.2b) 

Assuming that fields E(f,t) and H(f,t ) are varying harmonically with time, i.e. 
proportional to the oscillating term e lut , we can take the time derivative and next 
apply the curl operator to both equations (2.1a) and (2.1b). Next plugging VxH 
term into second equation (2.1b) and V x E term into first equation (2.1a) we split 
Maxwell’s equations into two subsystems, separating E and H fields: 


and 


V x V x E(r) = e(r)kgE(r ), 

V • (e{f)E{f) j = 0, 


V x V x H(r) 
e(r) 

V • H(r) 


k 2 0 H('r), 

0, 


(2.3) 


(2.4) 


where k Q = - = , and u> = J- — 

Each system of equations (2.3) and (2.4) contains the full description of electro¬ 
magnetic held. We can either use system (2.3) to find E held and next express held 
H in terms of V x E by using (2.1b), or alternatively, the held E[ can be obtained 
from equation (2.4) and held E can be found in terms of V x H with help of (2.1a). 

Here we prefer to use system (2.3) over (2.4) to avoid complications caused by the 
V x x H(f) term. 
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Now considering the vector identity [41]: 


VxVx£ = V(V£) - X7 2 E, 


(2.5) 


and expansion of the Gauss equation (2.2b): 


we obtain the following identity: 


V ( eE) = eVE + E ■ Ve = 0, 


Ve 


V x V x if = —V [ E ■ — ) — V 2 E, 


Thus system (2.3) can be written as: 


V‘E + tkiE = -V ( E ■ ^ 


( 2 . 6 ) 


(2.7) 


( 2 , 8 ) 


The equation (2.8) is all that is needed to describe electromagnetic waves in inho¬ 
mogeneous nonmagnetic media. Once the electric held E is found the corresponding 
magnetic held can be expressed in terms of V x E as follows from (2.1b) equation. 


2.2.1 Weakly guided approximation 

Let us assume that the rate at which function e(r) is changing with the change 
of r coordinate is insignificant, i.e. Ve(r) is smaller than e(f). Then the right hand 
side of equation (2.8) can be neglected [42]: 

V 2 E + e(f)k 2 0 E = 0 (2.9) 


This is the so-called weakly guided approximation. The equation (2.9) provides full 
description of electromagnetic waves in the case of small Ve(F). 

Let us write equation (2.9) in cylindrical coordinates: 


/v 2 


+ ekl 


& 


0 


0 

0 


V 2 + ek 2 0 ) 


\ (e p \ 


E ( 

\E Z J 


= 0 . 


( 2 . 10 ) 
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Here V 2, 0 is the Laplacian for a scalar function ifj in cylindrical coordinates: 

vy = + -a,+j 2 di +ad v, (2.H) 

and V 2 A is the vector Laplacian in cylindrical coordinates: 

V 2 H = ^ ~d p (pd p ) + + <9 2 ^ (a p p + A^cj) + , (2.12) 

The equation (2.10) can be obtained from (2.12) by taking into account that unit 
vectors themselves are functions of coordinates, and noting that in cylindrical coor¬ 
dinates there are two nonzero derivatives of the unit vectors: d^p = (j) and d^cj) = —p, 
the remaining derivatives are zero. 

We note that equation (2.10) decouples into two independent equations, 
the vectorial equation: 



and the scalar equation: 


X/ 2 E z + ek 2 E z — 0. 


0, 


(2.13) 


(2.14) 


Here we are interested in circularly symmetric optical waveguides, assuming that 
the waveguide structure is uniform along the z coordinate, and searching for a periodic 
solution in 0, i.e. u(p,(f>,z ) ~ u{p)e^ z e ^ m ^, we can can replace d z with j/3 and d$ 
with jrn. Thus the equations (2.13) and (2.14) can be written in form of the vector 
eigenvalue problem: 


K + H - A 1 + ek ° ) (A) 

V A + ekl) \eJ 

and the scalar eigenvalue problem: 

U f + '-d„ -A- A) E, = PE Z . 



(2.15) 


(2.16) 


Here A = (3 2 are unknown eigenvalues and /3 is the propagation constant. 
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The equation (2.16) is well known and is often considered for the weakly guided 
approximation case while the equation (2.15) is rarely used, see for example [43, 44, 
45], 

Once any of these two equations are solved, the remaining components of the 
electric field E can be found from Gauss’s law (2.2b), and next the magnetic field H 
can expressed in terms of electric field with use of equation (2.2b). 


2.2.2 The exact solution for cylindrical waveguides 

The exact solution can be obtained if term Ve is no longer neglected. The elec¬ 
tromagnetic field in such a case is described by equation (2.8): 


w - »/."/■: - -v (/■: Am 


(2.17) 


Let us assume that dielectric permittivity varies only along the the radial direction 
e = e(p), then 

p Ve _ d p e(p) 1 <9 0 e(p) d z e(p) _ 

e e(p) > e{p) e(p) 


= E, 


<V(d) 


" <P) 

= E p - (lne(p))', 


(2.18) 


~^ drl 


and 

/ (In e)"E p + (In e)'d p E^ 

V r' t) = V (^( lne )') = (In eyfaEp 

\ (In e)'d z E p y 

Now equation (2.17) can be rewritten in cylindrical coordinates: 
/V 2 — -4 + ekl + (In e)" + (In e)'d p 




^ + (lne)'i^ 


V 2 -4 + efc 2 


0 

0 


V 


(In e)'d z 


V 2 + efc 2 y 


\ (e p \ 

E<j, 

\E Z J 


(2.19) 


= 0. (2.20) 


The first two equations can again be separated, as they do not depend on the E. 
component. 
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Replacing d z with j/3 and d$ with jrn as previously, we obtain a slightly different 
from of equation (2.15). The vectorial eigenvalue problem (2.21), now depends on 
the rate of change of dielectric permittivity e'(p): 


Ul+ld p +(lneyd p -^+ekl + (lne)" \ fE p \ fE p \ 

V ^ + ?( lne )' d l + K- ^ \ E J \ E J 

( 2 . 21 ) 

Here we note again that terms e'(p) were previously ignored in the weakly guided 
approximation (2.15). In the following section we compare the results following from 
the weakly guided approximation (2.15) and from the exact formulation (2.21). 

Finally, once the fields E p and are known from equation (2.21), the remain¬ 
ing E z component can be found from the Gauss’s equation (2.2b) written in cylindrical 
coordinates: 

E, = --J~(^E p ) + —eX ( 2 . 22 ) 

j/3 \epdp p J 

and the magnetic held H can now be expressed in terms of electric held E = (E p , E^, E z ) 
with help of equation (2.2b): 

H = j—V x E. 

Ojpo 


(2.23) 
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2.2.3 TE and TM modes in slab waveguides 


/v 2 + ekl 

0 0 \ 

( Ex ) 


^(ln e)"E x + (In e)'d x E x ^ 

0 

V 2 + e/c 2 0 

Ey 

= - 

(In e)'d y E x 

V 0 

0 V 2 + eh 2 J 

\E Z ) 


V (In e)'d z E x j 


Here we show that our approach is valid, by deriving the well known equation for 
slab waveguides. Again, let us start with equation (2.8): 

V 2 i? + ek 2 E = —V ■ (2.24) 

Assuming that the slab waveguide is uniform along the y and z axis, the dielectric permit¬ 
tivity can only be function of the x coordinate e = e(x), thus equation (2.24) can be written 
in Cartesian coordinate system as follows: 


(2.25) 


Considering the waveguide symmetries, uniformity along the y and z axis, and as¬ 
suming that the wave is propagating along the z axis, we can can replace d z with j/3 
and d y with 0. The scalar Laplacian V 2, 0 can now be written in the following form: 

vV = (% + % + %)1> = 

= (dl - /3 2 ) Vb (2.26) 

and equation (2.25) can be rewritten in the the form of vectorial eigenvalue problem: 


(2.27) 


The first 2 equations are not coupled and can be written separately: 

(4 + (lneR + (hie)" + ckl) E x = /3 2 E x (2.28) 

(4 + ekl) E, = pE t (2.29) 

The equation (2.29) is well known equation for the TE modes in slab waveguide [46, 
47], 

The equation (2.28) can be recognized as equation for TM modes if we rewrite it 
in a slightly different form: 


^d 2 x + (In e)'d x + ek 2 + (In e)" 

0 

0 N 

( E x\ 


(e x \ 

0 

d 2 x + ek 2 0 

0 

Ey 

II 

ta 

to 

Ey 

V (In e)'P 

0 

d 2 x + ek 2 oJ 

\ E J 


U/ 
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where we have considered that 



u" + (lne)V + (In e)"u (2-31) 


We note that in the case of slab waveguide the system of equations (2.27) for 
the transverse field components E x and E y decouples into separate equations, inde¬ 
pendently describing TE and TM modes, unlike in the case of cylindrical waveguide 
where the system of equations (2.21) is coupled. 

We also observe that in the case of weakly guided approximation the equations 
for TM and TE mode (2.30,2.29) become identical, as 


ed x -d x u(x ) ~ d 2 x u(x ), 


(2.32) 


e 


thus the degeneracy occurs, i.e. the two different set of eigenfunctions eE x (x) and 
E y (x) corresponding to the same eigenvalues. 
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2.3 The numerical method 

In this section we describe the numerical procedure we implemented to obtain 
a vectorial mode solution and propagating constants of a dielectric waveguide with 
circular symmetry. The waveguide can be of an arbitrary refractive index profile, i. e. 
it can consist of any number of radially stratified layers made of various materials 
including absorbing material such as metals with the dominant imaginary part in the 
refractive index. In the following sections these modes, in other words eigenfunctions, 
are used as a basis for a more general problem of circularly symmetric dielectric 
waveguide with a small perturbation along the z axis. 

2.3.1 The scalar modes 

Let us start our analysis with the weakly guided approximation and the scalar 
Helmholtz equation (2.16): 

(d 2 p + - p d p - y + ek 2 )j E z (p) = f3 2 E z (p). (2.33) 

Rewriting the above equation in a slightly different notation we get: 

--dp(pd p ) + U m (p) 

where 

9 

771 

U m (p) = k 2 0 n 2 0 (p) - ?L. (2.35) 

P 

The eigenfunctions are denoted as R™(p) and eigenvalues as known as propaga¬ 
tion constants in waveguide theory. It should be noted that for each index m a set 
of different solutions will be obtained, keeping this in mind, let us omit the index m 
from the notation for simplicity. 

First, let us build a matrix representation of equation (2.34). This can be achieved 
by introducing an uniform grid (pi, p2, ..., Pn) and representing the unknown function 
R(p) as a vector with components Rj = R{pj)- Next, applying the finite difference 


R*(P) = m 2 K(p), (2-34) 
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method we can rewrite equation (2.33) in the following form: 



(2.36) 


with the central difference approximation: 


cijRj—i A bjRj A CjRj-\.\ — 3 Rj , 


(2.37) 


where 


1 11 
h 2 pj 2 h ’ 



1 11 
h 2 pj 2h 


(2.38) 


Here n 3 = n(pj) is the refractive index prohle of a given structure. The refractive 
index can be a complex number to take into account absorbing materials. 

Thus equation (2.36) can be written in the matrix form: 


[L\R = P 2 R, 


(2.39) 


where [L\ is the sparse 3-diagonal matrix with bj on the main diagonal, a 3 on the 
lower sub-diagonal and c 3 on the upper sub-diagonal. 

Considering that equation (2.37) is the second order differential equation, the two 
proper boundary conditions should be imposed: 

1. The condition at infinity p —> oo. 

This condition can be imposed by introducing an extra point Rn+ i, outside the 
computational domain. Here we are interested in guided modes, thus we can 
assume that there is no held at infinity and simply impose the zero boundary 
condition. By choosing a sufficiently large computational window we ensure 
that the held value is negligible at the computational boundary. 


Rn+i — 0 . 


( 2 . 40 ) 
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Plugging this condition into (2.37) we have 

ajyR.N-1 + biyRjy + Cat • 0 = /3 2 Rn- (2-41) 

We note that this condition is already imposed by the initial equation (2.37), 
thus there is no need to modify this equation. 


2. The condition at p = 0. 


(a) Assuming that m — 1, ±2, ±3,... we should require 

R(p = 0) = 0. (2.42) 


By introducing an imaginary point f? 0 outside the computational domain 
on the left, we can write: 

cl\Rq T b\R\ T C\R 2 — /3 2 R\, (2.43) 

considering (2.42) we get: 

ai • 0 + b\R,\ + Cl R 2 = f3 2 Ri, (2.44) 


Again, we see that this condition is already imposed be the default form 
of the equation (2.37). 


(b) 


Assuming that m = 0 we get a maximum in the field distribution at the 
center of the waveguide, and thus should require a zero derivative of the 
field: 


dR (P ) | _ n 

dp lp=pl ’ 

using the central difference the above relation can be written as: 


(2.45) 


Ro = R-2 = 0 . (2.46) 

Plugging the above condition into (2.43) we get: 

b\R\ + (ci + ai)R 2 = P"Rij (2.47) 

Therefore the matrix [L] has to be modified: 

L 12 = L 12 + ai (2.48) 


This is the only modification of (2.37) which is required. 




Chapter 2. A Full-vector complex mode solver 


22 


Finally considering modification (2.48) the problem (2.37) can be stated in the 
matrix form: 

[L\R = P 2 R. (2.49) 

Here we note that the zero boundary condition can be imposed in both cases if the 
problem is solved on the internal p G [-R, R] instead of p G [0, R], 

The problem (2.49) is a well known numerical eigenvalue problem, which can be 
solved with variety of numerical techniques. Because we are dealing with a 3-diagonal 
sparse matrix, very efficient iterative methods for large sparse linear systems can be 
implemented. We used MATLAB software, which allowed us to consider up to one 
million elements. Even for a very large matrix of 10 6 x 10 6 elements, the required 
eigenvalues and corresponding eigenvectors, in the range of interest, were found in 
less than 1 minute. The matrix does not actually includes all of the 10 12 numbers, 
but rather 3 • 10 6 numbers written in terms of three separate vectors corresponding 
to the three matrix diagonals. The computational algorithm is based on the iterative 
routine converging from a random guess [48] towards the exact eigenvectors and 
corresponding eigenvalues in the given range. If, instead of a sparse matrix, a dense 
matrix is used, the computation might take a prohibitive amount of time for the 
particular problem presented in this work (for example, it would take a couple of 
days to compute eigenvalues of a 5000 x 5000 matrix in Matlab with the default 
computational routine). 

After the numerical eigenvalue problem is solved, a set of eigenvalues /3k and 
corresponding eigenvectors Rk are accessible. It is convenient to use a potential 
barrier analogy with quantum mechanics. The potential energy function is defined 
by equation (2.35). Basically, the square of refractive index profile n 2 (p) plays a role 
similar to the role of a potential energy barrier in quantum mechanics. The eigenvalue 
can be though of as “energy states”, and should be expressed in the same units as the 
potential barrier, or refractive index. The eigenfunctions Rk{p) represent the electric 
held profile, that is E z component if the scalar equation 2.33 is considered, or E p and 
E'</) components for the vectorial case, as will be discussed later. 

Here we introduce the so-called effective refractive indices i.e. N e ff^ = fcoAc, and 
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plot the normalized eigenfunctions Rk{p) centering them at the corresponding 7V e //^ 
with analogy to quantum mechanics. The result is shown in Figure 2.1. 




Figure 2.1: The refractive index profile of SMF-28 fibre immersed in water. The 
eigenvalues and eigenfunctions are plotted for m = 0 and rn — 1 

As can be seen from Figure 2.1 that in the case when N e ff t k becomes smaller than 
the refractive index of the surrounding medium, the held energy is leaking outside the 
fibre boundaries, and hence the energy is no longer confined inside the fibre as shown 
in Figure 2.2. Such modes are called the leaky modes. The proposed method should 
be modified for the proper treatment of leaky modes. The zero boundary condition 
should be replaced with a numerical absorbing boundary condition, ensuring that 
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the incident travelling waves are absorbed and no energy is reflected back into the 
waveguide, and hence the formation of standing waves is prohibited. Such absorbing 
boundary layer is called the perfectly matched layer (PML). An efficient numerical 
technique for implementing the absorbing boundary condition in the FDTD method 
was developed by Mur in 1981 [49]. 



Figure 2.2: The fibre cross-section and corresponding refractive index profile with 
modes bounded inside the fibre. 


We can note the further similarities with quantum mechanics. The number m 
in the equation (2.34) plays a role similar to the role of orbital quantum number in 
quantum mechanics. The larger the angular momentum of a particle, the closer its 
wavefunction is to the periphery, which is also consistent with classical mechanics. 
Unless the critical value is exceeded, the particle stays confined by central potential, 
although some fraction of its energy leaks outside the potential barrier boundaries. 
The described effect can be observed in Figure 2.3. The potential function U m (p ) is 
defined by the equation (2.35) and is shown in red color. We can note in Figure 2.3 
how significant is the role of angular momentum potential U m (p ) = 
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m = 40 



Figure 2.3: The radial eigenfunctions R™(p) (shown in blue color) and the potential 
well function U m (p ) (shown in red color) for m = 40. 

Finally, the complete basis functions ^?(p, </>) can be constructed by considering 
the angular dependence as well: 

W(p,<P) = K(py mt - (2.50) 

A particular solution can be constructed by choosing an orbital number m, thus 
defining the potential function U™(p), and picking a particular radial eigenfunction 
R";'(p) for the chosen mode family m. An example of such particular solution is shown 
in Figure 2.4. 


2.3.2 The vectorial modes 


Let us first consider the vectorial equation (2.15) for the weakly guided approxi¬ 
mation: 


'dl + i d p -^ + ekl 


P ' p* 

j2rn 


j2m 

P 2 


dl + -d 0 — ^ + ek 2 n , 


E p\ _ p2 ( E P 
Erh / \ Ea. 


(2.51) 


p* "P 1 p^P p 2 1 W W </V 

Again, we can implement the finite difference method to replace the system of two 
second order coupled ordinary differential equations (2.51) with a system of algebraic 
equations: 

[M\E r = (3 2 E t . 


(2.52) 
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Before constructing the finite difference implementation we note that matrix in (2.51) 
contains complex numbers. Although the complex numbers does not impose any 
limitation on the numerical method, the computational speed can be increased if the 
matrix is real. The memory requirement is also reduced by a factor of two for the real 
numbers. Multiplying the first equation in (2.51) by imaginary unit j and pulling it 
under the E p component we get: 

(2.53) 

<e f + \d f -^ + ekl)\Ej \eJ 

The operator matrix now is converted to a pure real form, unless the material permit¬ 
tivity e is complex. The same equation was previously obtained elsewhere [45, 43]. 

The finite difference matrix [. M] can be constructed in the following way: 

[m\e t = P 2 e t , 


'd 2 p +ld p 


m 2 +1 


+ eE 


2 m 


(2.54) 
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where 


[M] = 

L a \ 

( J, 

V« l ) 


2m 

Cii — 

„2 1 


Pi 

L = 

{d?+-m 


m 2 + 1 
P 2 


+ eki 


(2.55) 


here [ D ] 2 and [D] are the second order and the first order finite difference matrices, 
respectively, assembled with the help of the central difference approximation in a 
fashion similar to (2.37) and (2.38). We note that operator L becomes identical to 
the scalar operator (2.36) if m 2 + 1 is replaced by m 2 . 

The structure of matrix [M\ is shown in Figure 2.5, where the number of nodes 
was limited to N = 10 for the convenience of representation, however for a practical 
application at least several thousands of elements should be considered. 

In can clearly be seen that the matrix is sparse with only five main diagonals. 
The right most and left most diagonals represent term ^ which is responsible for the 
coupling between E p and E$ field components. In the case when m — 0 the matrix 
becomes a three diagonal, with uncoupled E p and E$ components. 

In the similar way we can approach the exact equation (2.21): 


K + \d p + (In e)'d p - + ekl + (In e)" f? \ (jE p \ = (jE p \ 

V ^ + 7M dl + \d p -^ + ekl) \eJ \eJ' 

(2.56) 

here we pulled the imaginary unit j out of the matrix. 

The system of algebraic equation takes the following form: 




Ex — (3 2 E t , 


(2.57) 


where a and L are defined previously in (2.55), and 

bi = ajH-(me) , 


P 

L 2 = L + (lne)'[D} + {\ne)". 


(2.58) 






Chapter 2. A Full-vector complex mode solver 


28 


* * 

* * ♦ 

* * * 


* * * 

♦ * ♦ 


♦ * 

* * * 


♦ ♦ 

♦ ♦ ♦ 

* * * 

* * * 
♦ * 


* * ♦ 
♦ ♦ 


10 


15 


20 


Figure 2.5: The structure of the sparse matrix [M] for N — 10 and m ^ 0. 


The matrix for the exact case has the same structure as shown in Figure 2.5. 

Alternatively, the exact system of equations (2.21) can be rewritten in the Sturm- 
Liouville form: 

i (? + ?(>" A) pf„ (Ji) + (-f + <) 

(2.59) 


jepE P \ = 02 (:') ( P E i\ 

, peJ \ P eJ ’ 


= u" + -u — —u, 

= u" + -u' + (lne)V —+ (lne)"u. (2.60) 

We discuss this subject in more detail in the following section. Next, the numerical 
finite difference method explicitly designed for Sturm-Liouville problems can be ap¬ 
plied. The main idea of this approach is to derive a central difference equation for 
the T P (piP)^) operator. 

In the next section we compare results obtained by numerically solving the scalar 
(2.14) and vector (2.13) equations for the weakly guided approximation with the exact 
solution obtained from (2.21). 


where we assumed 
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2.4 Discussion 

In this section we verify onr results and discuss properties of the inodes. First let 
us start by verifying the developed numerical method. The analytical solution can 
be easily obtained for slab waveguides, as shown for example in [46, 50]. 

The solution for TE and TM modes guided in a thin symmetric glass slide im¬ 
mersed in water as well as the dispersion curves were presented in [46]. 



Figure 2.6: The dispersion curves of a symmetric glass slide immersed in water. 
The waveguide and the medium refractive indices are rtwG = 1-45 and riMed. = 1.33, 
respectively. 


Defining the normalized frequency V and normalized propagation constant b as: 


V 

b 


Rk 0 

N 2 

ly eff 


Tin 


R 

A’ 


n \ — 77 .| 


N 2 

7 V e//’ 


(2.61) 


and solving equations (2.28, 2.29) numerically using the proposed method we plot the 
dispersion curves for TE and TM modes, as shown in Figure 2.6. The identical graph 
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was obtained in [46]. Here R is the waveguide width (or diameter), ri\ and n 2 are the 
refractive indices of the waveguide and the surrounding medium, respectively. In the 
similar manner we can verify the results for cylindrical waveguides, see for example 
results obtained elsewhere [17, 18, 44, 45, 46, 50]. 

It was shown in the previous section that the Maxwell equations can be decoupled 
into the scalar (2.14) and vector (2.13) equations if the weakly guided approximation 
applies. Solving these equations independently we obtain dispersion curves as shown 
in Figure 2.7 and Figure 2.8. 

The eigenvalues for the scalar and vectorial problems are identical. Indeed, this 
result can be expected as both equations were obtained from the same initial system 
of Maxwell’s equations, and are in a way simply a different representation of the same 
physical system. 

We note that solutions to the vectorial equation are degenerate, except for the 
modes with the zero angular part, i.e. m = 0. The equations for the E p and E$ field 
components become uncoupled at m — 0, thus E p component can take arbitrary val¬ 
ues independent on the value of E^ component. In other words the modes belonging 
to the m — 0 family can have an arbitrary polarization, with no restriction imposed 
by the cylindrical symmetry of the problem. In such cases modes can be represented 
either by a superposition of two orthogonal linearly polarized waves or as a super¬ 
position of left and right circularly polarized waves. In the literature these modes 
are known as linearly polarized (LP) modes [17, 43], although the term circularly 
polarized (CP) modes would also be appropriate [43]. 

The mode profile at m — 0 is shown in Figure 2.9. Indeed, it can be noted that 


two modes 


E P 

0 


and 



have an identical propagation constant, thus can be su¬ 


perposed to represent an arbitrary polarization, including the circularly polarization. 
Hence the widely used term the linearly polarized modes are in a way misleading. A 
more detailed discussion on this subject can be found in [43]. 
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Figure 2.7: The dispersion curves of the scalar (2.14) modes, here riwa = 1-45 and 
nMed. = 1.33. 



Figure 2.8: The dispersion curves of the vector modes (2.13). The waveguide and 
medium refractive indices are nwc = 1-45 and riMed. = 1-33, respectively. 
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a) (i,i) 



p. pm 


b) 




IE/ 
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Figure 2.9: The degeneracy of modes: the scalar (1,1) mode and two vectorial (0,1) 
and (2,1) modes. The waveguide and medium refractive indices are uwg = 1-45 and 
nMed. = 1-33, respectively. The potential barrier is depicted with green color 
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In addition we note that the first mode in a cylindrical waveguide, or the core 
mode in a single mode fibre, is obtained by assuming that rn — 0 in the scalar case or 
m — 1 in the vectors case, as can be seen from Figures 2.7 and 2.8. In the vectorial 
case the second mode is observed at m = 0 , with which the so-called LP modes are 
associated. Therefore the LP modes can propagate only in the cladding of a single 
mode fibre. 

In the case when the weakly guided approximation can be applied, each mode 
can be obtained by solving either the scalar equation (2.14) or vector equation (2.13). 
The vectorial modes are mostly degenerate. It is interesting to compare radial profiles 
of such identical modes. For example, in Figure 2.9 the radial profile of (1,1) scalar 
mode and (0,1), (2,1) vectorial modes are shown. All three modes have an iden¬ 
tical propagation constant. We note that the LP modes (0,1), with the radial and 
angular components decoupled, coincide with the (2,1) mode, with coupled E p and 
E^ components. It is interesting to note that all the modes have the identical radial 
profile, even though the scalar mode was obtained by solving the scalar equation at 
m = 0 for the E z component, where the vectorial modes were solved for E^ and E^ 

(E p \ 

at m = 1 and m = 2. Each mode is normed to unity, hence the difference in 


. Ea 


amplitudes between (0,1) and (2,1) modes. The (0,1) mode has only single nonzero 


component: 


E n 


or 


0 

. Ea 


, whereas in the mode (2,1) the both components E p 


and E^ are nonzero. 

Let us consider the vectorial modes with m/0. The held components E p and 
E<p are always coupled in this instance, thus one can not be chosen independently of 
another. In the case of weakly guided approximation the modes are mostly degenerate. 
For example, the radial profile of (2,2) scalar mode and (1,4), (3,3) vectorial modes 
are shown in Figure 2.10. Again, we note the similarity in the radial profile of these 
modes. 

If an exact solution is considered the degeneracy is removed. The modes with a 
different m number start to behave differently, diverging apart from each other with 
the increase in refractive index difference between the waveguide and the medium, as 
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( 2 , 2 ) 


' E zl 


2 



p, p m 



Figure 2.10: The degenerate modes: the scalar (2, 2) mode and two vectorial (1,4) and 
(3,3) modes. Here the waveguide with n -2 = 1.45 is immersed in water riMed. = 1-33. 


shown in Figure 2.11. The modes (1,2) and (3,1), (1,4) and (3,3), (1,6) and (3,5) 
no longer coincide, as in the case of weakly guided approximation. In the same graph 
the single modes (1,1), (1,3) and (1,5), without a pair mode, are also plotted. In 
general, each mode obtained for a small refractive index difference, splits into two 
modes, each with a distinct dispersion curve. 

From Figure 2.11 we note that not only the dispersions curves of different mode 
families m diverge apart from each other, but even the modes belonging to the same 
family start to behave differently, for example the (1, 3), (1, 5) and (1, 7) modes. This 
effect can be studied in more detail by plotting the dispersions curves corresponding 
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Figure 2.11: The split in dispersion curves between m = 1 and m — 3 modes. 

The waveguide and the medium refractive indices are nwc = 3 and iiMed = 1-33, 
respectively. 


to different refractive index differences An = nwc — nMed between the waveguide and 
the medium, as shown in Figure 2.12 for the rn — 2 family. 

We note that different modes inside the same mode family m, behave differently. 
For example, the modes (2,2) and (2,4) have barely changed in location on the 
dispersion plane, whereas modes (2,1), (2,4) and (2, 5) are highly divergent, especially 
at the cutoff region. We discuss this effect for the rn — 0 case in more detail later. 

We conclude our discussion here by pointing that if m ^ 0, and the refractive index 
difference between the waveguide and the medium is significant, so that the weakly 
guided approximation can no longer be applied, the majority of modes splits into two 
separate modes, thus the degeneracy is removed. Once this property is understood 
we can move on to a more complicated case of m = 0. 

If the vectorial equation is solved for the case of m = 0 the degeneracy observed 
for a small refractive index difference is removed as was discussed above. But now 
not only the modes with different m numbers are distinctly seen as separated (for 
example, as was shown in Figure 2.11 ), but additionally the modes inside the same 
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( 2 , 1 ) 

( 2 , 2 ) 

( 2 , 3 ) 


( 2 . 4 ) 

( 2 . 5 ) 


Figure 2.12: The dispersion curves, obtained by solving the exact vectorial equation, 
for various refractive index ratios An = n W c ~ n Med. at m = 2. The waveguide is 
immersed in water UMed. = 1-33. 


family m = 0 are split into two subfamilies. Thus a single mode is split into three 
distinct modes, as shown in Figure 2.13. 

The split between modes belonging to different families m — 0 and rn — 2 was 
discussed above. Now, we focus our attention to the extra split within the same 
family of modes at m = 0. To understand this split we shall go back to the vectorial 
equation (2.21): 

(dl + )d p + (In e)'d p -^±i + ekl + (In e)" \ f E p \ = (E p \ 

V 1E + dl + H p -^ + ekl) \eJ P \eJ ’ 

(2.62) 

which splits into two uncoupled separated differential equations at m = 0: 

{^d 2 p + —d p + (In e)'d p + ek~ + (lne) ,/ ^ E p = fd 2 E p , (2.63) 

(yd 2 p + —d p + E<f, = f5 2 E 


and 


(2.64) 
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Figure 2.13: The split in dispersion curves between m — 0 and m — 2 mode families, 
and split between the modes inside the same family at m = 0. The waveguide and 
the medium refractive indices are nwc — 3 and nued. = 1-33, respectively. 


In the weakly guided approximation case the term (In e)' can be neglected, thus 
equation (2.63) becomes identical to equation (2.64), hence the same propagation 
cases are obtained from both equations, and hence the corresponding modes are 
degenerate. 

However, if the refractive index difference between the waveguide and the medium 
is significant, the term (lne)' can no longer be neglected, and thus different propa¬ 
gation constants would result from the solution to equations (2.63) and (2.64). The 
initial degeneracy would be removed. 

The result for two modes at m = 0 is shown in Figure 2.14 for a significant re¬ 


fractive index difference. The modes 


E P 

0 


and 



are obtained by solving equa¬ 


tion (2.63) and (2.64), respectively. The split, between the corresponding propagation 
constants is clearly seen. Thus the modes no longer coincide and hence can no longer 
be used to represent an arbitrary polarization, such as linear polarization, therefore 
the typically used term “LP modes” does not completely reflect the physical proper- 
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ties of such modes. 


m = 0 
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Figure 2.14: The split between the eigenvalues of pure radially polarized E p and pure 
angular polarized modes of the nn — 0 family. The potential barrier is depicted 
with green color, the core refractive index n 2 = 3 is surrounded by the cladding with 
TT-i = 1.33. 


The vector E of the 


E P 

0 


mode is orientated radially, i.e. is aligned transver- 


sally with respect to the interface along the p vector, whereas the vector E of the 
. 0 , 

mode | I is aligned tangentially to the interface, along the <f> vector. The remaining 


,E, 


<t>; 


E z component gives a small contribution to both radial and tangential components. 

The equations (2.63) and (2.64) can be rewritten in a more elegant Sturm-Liouville 
form: 


(p-T- (-J-) + ek l) v = p 2v -> (2.65a) 

V dp \pdpj ) 

( ep ~r (—~r) + ek °) u = ^ u ' (2.65b) 

V dp \pedpj ) 


Here v = pE$ and u = jepE p . The equations were derived with the help of (2.60). 
We note the similarity between equations (2.65b, 2.65a) and the eigenvalue equa- 
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tions for the TE 2.29 and TM 2.30 modes in the case of slab waveguides: 


(J? + ek ^j v = P 2v (2.66a) 

(4 ( 4 ) +ek ) “=^ (2 - 66b) 

Here v = E y and u = eE x . Because of this similarity the solutions obtained 
form (2.65a) and (2.65b) equations sometimes called TE and TM modes, respec¬ 
tively. The role of E x and E y components in a slab waveguide is played by E p and 
E^ components in the case of a cylindrical waveguide. 

Both equations (2.66b) for the TM modes in the case of a slab waveguide and 
the equation (2.65b) for the TM-like modes in the cylindrical waveguide case contain 
terms and ep respectively, are sensitive to the rate of change of 

permittivity. Therefore, we conclude that the TM and TM-like modes should be 
particular sensitive to changes in refractive index. 



Figure 2.15: The dispersion curves, obtained by solving the exact vectorial equation, 
for various refractive index ratios An = 71wg ~ n Med. at m — 0. The waveguide is 
immersed in water n = 1.33. 
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The dispersion curves for various index contrasts are shown in Figure 2.15. Indeed, 
it can be clearly seen that the TM-likc modes for the E p component are seen shifting 
significantly with the increase in index contrast, whereas the TE-likc modes for the E$ 
component are almost immovable. The operator ep^ sensitive to changes 

in the refractive index profile, affecting the E p solution. On the other hand, the 
rate of change in the refractive index is not incorporated in the E^ solution, hence 

_/y2 _ n 2 

if the corresponding propagation constant is plotted along the b = ^/_ n 2 2 axis the 
corresponding dispersion curve should stay intact. 
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2.5 The orthogonality of the basis functions 

In the next Chapter we are going to look for a solution to the more general problem 
of a waveguide with a small perturbation along the z axis. The modes obtained for 
the unperturbed case might be used as a basis function in terms of which the general 
solution can be expressed. The calculations can be significantly simplified if the basis 
functions are orthogonal. In this section we derive the orthogonality relation in terms 
of a weighted product, ensuring the orthogonality of modes. 

To our knowledge the orthogonality relation for the vectorial modes of the form 
- (E p \ 

E = in cylindrical waveguides with an arbitrary radial refractive index profile 

\eJ 

is not available in literature. 

2.5.1 The orthogonality relation for the scalar modes 

Let us start by checking the orthogonality property of scalar modes. The equa¬ 
tion (2.16) 


d 2 Id m 2 , 9 i „ 

d? + - p Tp-y +Ez = p E > 


can be rewritten in the self-adjoint or Sturm-Liouville form [51, 52]: 

e(p) = 0, 


r ( p ^7 ~) + q{p) + Md) 

dp \ dp 


that is 


where 


d f d \ / , 9 rrv 

Q 


Tp{%) + { ek ‘° p -T ] ~ f)2f ’ eW = 0 ’ 


(2.67) 


( 2 . 68 ) 


(2.69) 


P(p) 

q(p) 
w(p ) 


= P 


t.2 m 

tkoP - 

p 


= p 


(2.70) 


It can be proven that the eigenfunctions ek(p) obtained by solving equation (2.68), 
written in the Sturm-Liouville form, are orthogonal with respect to the weighted 
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function w(p), i.e.: 

<e k \w\ej>= w(p)e k (p)ej(p)dp = 8 kj , (2.71) 

J a 

where e k (p) and X k are the eigenfunctions and eigenvalues of (2.69), respectively. 

Here we provide a short proof of the fact that the scalar modes obtained from 
equation (2.68) are indeed orthogonal, and later we will build our original prove for 
the vectorial case upon it. Detailed discussions of the scalar case can be found in [51] 
and [52], although a more general proof for the vectorial case was not provided. 

We start with the eigenvalue problem written in the Sturm-Liouville form (2.68): 

Le k (p) = X k w{p)e k {p). (2.72) 

where 

r d f d\ 

- Tp Ydp) + q - 

Considering 

Lu = pwu , 

Lv = Xwv, 

We can construct the following expression: 


(2.73) 


(2.74) 


vLu — uLv — (ji — X)wuv. (2-75) 

Inserting the operator L and using the Lagrange’s Identity [51, 52] we obtain 

dp{pud p v — pvdpu) — (/i — A )wuv. (2.76) 

After integration over the interval p G [a, b] 


/ dp{pud p v - pvdpu) =/(//- A )w(p)u(p)v(p)dp, 

J a J a 

we come to the Green’s Identity [51, 52]: 

P{P) = ^ ~ ^ L w ^ u ^ v ^ dp - 


(2.77) 


(2.78) 
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The left hand side vanishes if we consider that the field decays outside the waveguide 
boundary and is approximately zero at the numerical boundaries. 


u(a) = v(a) = 0 , 
u(b) = v(b ) = 0 . 


Finally we obtain 

0 = {n ~ A) / w(p)u(p)v(p)dp. 

J a 

For non-degenerate modes the above relation can be true only if 

f b 

/ w(p)u(p)v(p)dp = 0 , for p 7 ^ A. 

J a 

Thus we have proven the orthogonality relation for the scalar case (2.69): 


< E Z JE Z>X >= / pE Ztli (p)E ZtX (p)dp = 0 , for p ^ A. 


(2.79) 


(2.80) 


(2.81) 


(2.82) 


J a 

In the next section, in a similar manner, we prove the orthogonality relation for the 
vector case. 


2.5.2 The orthogonality relation for the vectorial modes 


Let us assume that we have a system of equations: 


L\ 

b 



wu 

0 



(2.83) 


where operators L i and L 2 are given in the Sturm-Liouvillc form: 

Ll = T P { pM ^) +qM ’ 

^ 2 = ~r ( P 2 ^j~) + q 2 ^- ( 2 - 84 ^ 

dp \ dp J 

Following the derivation in the previous section let us start with two modes: 


[M]u = X[W]u, 
[M]v = p[W]v. 


(2.85) 
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Our goal here is to find the orthogonality relation between the two modes u and v. 
First, let us construct the scalar product: 

tf[M]u = \tf[W]u, (2.86) 

here ip = iy 1 )*■ Rewriting the above relation in coordinate form we obtain: 


(O Ob, 


L\ g 
h L 2 


u l 

02 , 


= \(v*i,v* 2 ) 


W 11 0 

0 w n 


u 1 

02 , 


or 


v\LiU\ + gv\u 2 + hv 2 u\ + V2L 2 u 2 = A v^[W]u. 

Now the expression 

v ! [M}u — u\M]v = A v ! [W]u — pu^[W]v 

can be rewritten in the form 

(ui u 2 - u 1 v 2 )(g - h) + (y\L\U\ + v 2 L 2 u 2 ) - (uiLy^ + u 2 L 2 v 2 ) = 

= (A — n) (w 11 u 1 v 1 + w 22 u 2 v 2 ), 

or rearranging the terms we get: 

(viLiUi - RiZ/iVi) + ( v 2 L 2 u 2 - u 2 L 2 v 2 ) = 

= (A - g)(w u uiVi + w 22 u 2 v 2 ) + (ui u 2 - u\v 2 )(h - g). 

Finally, applying the Lagrange’s Identity we obtain: 


/ du\ dv i 


+ 


, , du 2 dv 2 

Mp] ( - M2 — 


= (A - p) / (wn(p)ui(p)ui(p) + w 22 (p)u 2 (p)v 2 (p))dp + 


(MpWp) -«i(p)o(p))(%) -g(p)))dp, 


(2.87) 


( 2 . 88 ) 


(2.89) 


(2.90) 


(2.91) 


(2.92) 


Assuming as previously that at the computational boundaries all the field components 
vanish, we get the following expression 

(A — p) [ (wii(p)ui(p)ui(p) + w 22 (p)u 2 (p)v 2 (p))dp = 


= / (Mp)« 2 (p) -«i(p)o(p))(%) -g{p)))dp. 


(2.93) 
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For the two different modes, i.e. n ^ A, the expression 


< u\W\v >= / (tun(p)tti(p)ui(p) + w 22 (p)u 2 (p)v 2 (p))dp ^ 0 


(2.94) 


is non zero, thus we can not consider inodes u and v to be orthogonal in the gen¬ 
eral case. However, as we will see later, for a small azimuthal number m we can 
proximately consider modes to be orthogonal: 


< w|hF|?; >~ 0 , for p, ^ A. 


(2.95) 


Now, let us apply the orthogonality expression (2.94) to a particular case of modes 
propagating in cylindrical waveguide of an arbitrary profile (2.21): 


d 2 p + i d p + (In e)'d p -=**! + ekl + (In e)" - P 2 
i%L + im.a ne y 

P P X J 


_j 2 m 

P 2 


dl + l p d p -^ + ekl - p 2 


Ep \ = 0 (2.96) 

ty j 


First, let us rewrite equation (2.96) in the Sturm-Liouville form (2.83). Consider¬ 
ing that 

1/ vY „ 1 , 1 

-(to) = u H —u -TO, 


( —(reu)'^ = u" + to' + (In e)'u' - + (In e)"u (2.97) 

\re J r r z 


we can rewrite the matrix (6.24) equation in the following form: 


1 j 2 m 


(u P pt) + ^{-f + <-p 2 ) pt 

f (W) i 0T P P ) + \(rf + ekl - P 2 ) p\ \e, 




= 0 (2.98) 


Now we can pull ep and p under the E p and held components, respectively: 


L\ g 
h L 2 


f tpEp 
\P E <P 




ep 

0 A 

P J 


f e pE p 

\P E P. 


(2.99) 


here 


g(p) = 


1 2 jm 
P P 2 


7 / \ 1 2 jm jm 

Kp) = — (—2"+ lne - 

ep V P P 


( 2 . 100 ) 
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f d [ 1 d\ 1 [ m 2 2 

- TpXJeTp) + y + tk ° 


d f 1 d \ 


m 


i2 - Tp \j,tp) + ~p + ek °) (2 ' 101) 

Thus, the eigenvalue equation for cylindrical waveguide (2.96) can be written in 
the form of equation (2.83). Considering that 

, , , , . 1 (2jm ,, jm\ 12 jm 

h(p)-g(p = - -V + One — +--— = 
e p V p 2 p J p p 


pe v P 

a standard orthogonality relation follows form (2.93): 


jm ( 2 e _±l + (^ e yi\ = o for m = 0, (2. 

V P PJ 


102 ) 


< u\W\v >= / (w u (p)ui(p)vi(p) + w 2 2 {p)u 2 (p)v 2 (p))dp = 0 , for p ^ A. (2.103) 
J a 

The equation (2.103) approximately holds for small azimuthal numbers m, as we will 
see in the next section. The above expression can be rewritten in terms of vectorial 
f E n 


modes E = 




in cylindrical waveguides as follows: 


< u\W\v >= / (wii(p)ui(p)ui(p) + w 22 (p)u 2 (p)v 2 (p))dp = 
J a 


= [ ( — {peE p ^y(peE PtX ) + -{pE^)*(pE^ x )\dp (2.104) 

J a V P^ P J 

Hence, for the same mode family with a relatively small azimuthal number m : 

< E^Ex >~ [ P {t*E* p4l E P) x + E^x) dp = 0 , for p ± A. (2.105) 

J a 


The above equation is used extensively in the following sections. We note that the vec¬ 
torial orthogonality relation (2.105) resembles the scalar orthogonality relation (2.82). 
The same weighted function p is used in both cases, although the extra weighted func¬ 
tion e(p) is also used to weight the E p component. 


2.5.3 Numerical verification 

Solving the eigenvalue problem (2.16) or (2.21) numerically we obtain a set of basis 
eigenvectors (efc(p)}- We can verify the solution by checking whether the obtained set 
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of eigenvectors satisfies the orthogonality relations (2.82) for the scalar case (2.16), 
or (2.105) for the vector case (2.21). 

We note that for a different m number we get a different Sturm-Liouvillc problem, 
with a different set of eigenvalues {/d™} and eigenvectors {efc(p) m }. For example, for 
the scalar case (2.16) we have: 


dp(pdp) + 




Kp 


K(p) = o, 


(2.106) 


here index m denotes the family of modes obtained for a given angular momentum 
defined by the m number. 

According to the equation (2.82) the orthogonality relation is defined as 


noo 

< R£\R? >= / pR n k {p)R™{p)dp, (2.107) 

Jo 

here we considered the proper boundary conditions on [a, b] = [0, oo]. 

Thus equation (2.106) allows us to construct a complete set of basis functions, 
orthogonal in the sense of (2.107). As we noted previously, the orthogonality property 
will substantially simplify the analysis of a more general problem in the following 
chapter. 

Now let us construct the overlap matrix | R™ >< R^ \ = CfjJ 1 representing the over¬ 
lap integral between j'-th and k-th modes, taken from the m-th and n-th families. 


/»oo 

\R?Xl%\ = pK(p)RJ l (p)dp = Cjr = 0 if J^k. ( 2 . 108 ) 

Jo 

For convenience, let us norm each basis function to unity: 


p\RT(p)\ 2 dp = l. 


(2.109) 


We can view the normalization procedure as a construction of a new basis: 

KXp) 


RTXp) = 


sj 7o°° P\ R kP)\ 2 dp 


( 2 . 110 ) 


The resulting overlap matrix CJ™ is shown in Figure 2.16. As can be seen, all 
the modes belonging to the same family of modes, with the identical m number, are 
mutually orthogonal, regardless of the refractive index profile n(p). However, if the 
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overlap integral is computed between modes of different families, for instance R k (p) 
and Rj°(p), the result is non zero in most of the cases. The overlap matrix is dense. 
Indeed, each family of modes is a solution to a different Sturm-Liouville problem, 
uniquely defined by a different potential barrier (2.35). 


m=0 to m=0, or m=10 to m=10 



The mode number 


m=0 to m=10 

0.6 
04 

■ • 0.2 
■ ■ 0 
- 0.2 
- 0.4 

U - 0.6 

: 

5 10 15 20 25 30 35 40 

The mode number 



Figure 2.16: The overlap matrix C”™ 1 is computed for SMF-28 fibre immersed in 
water. The matrix is an identity for modes belonging to the same family of mode 
(qr = I for 77i = 0 or m = 10). However, if a different families of modes (m — 0 
and n = 10) are considered, the overlap matrix C is dense. 


A similar analysis can be conducted for the vectorial case. The orthogonality 
relation for a small azimuthal number m can be written in accordance with (2.105) 
as follows: 


qr = \R? >< K\ = / + 0 *f ( 2 . 111 ) 


The radial function can be normed for convenience: 

q m (p) 


H?(P) = 


d p 


( 2 . 112 ) 


The overlap matrix q™:” is shown in Figure 2.17. 

We note that for the vectorial case the number of modes is approximately doubled 
in comparison with the scalar case due to the mode splitting, as was discussed in 
Section 2.4. 
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m = 0tom = 0 m = 3 to m = 3 




The mode number The mode number 
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m = 1 to m = 0 


m = 1 to m = 2 
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Figure 2.17: The overlap matrix for the vectorial case is computed for SMF-28 fibre 
immersed in water. The overlap matrix C^, m , for the modes belonging to the same 
family m, is an identity matrix for m = 0 or close to the identity for small azimuthal 
numbers. For different families of modes, the overlap matrix is dense. 


As can be seen from Figure 2.17, the orthogonality relation holds almost exactly 
for small azimuthal numbers m. In the next section we will see that we need modes 
only from a few first families with a small m number. In addition, we will show that 
due to the phase matching condition, the overlap within a family of modes can be 
neglected, except for the m = 1 family. 

We conclude that radial components of the basis functions are not mutually or- 
















Chapter 2. A Full-vector complex mode solver 


50 


thogonal for different families of modes, i. e. families with different azimuthal numbers, 
in both scalar and vectorial cases. However, within a particular family m the radial 
components can be considered to be mutually orthogonal. In the following section 
we will see that although radial components belonging to different families of modes 
may not necessarily be mutually orthogonal, the complete basis functions are always 
orthogonal due to the orthogonality of angular components. 

2.6 Conclusion 

We conclude this section by stating that the scalar model yields the same result 
as the vector model if the refractive index contrast is small and hence the weakly 
guided approximation can be applied. However, if the index contrast is high each 
scalar mode splits, in general, into two separate modes. 

In the vectorial case the first mode, with the highest effective refractive index, 
occurs at m = 1 and is single. The remaining modes come in pairs (for exam¬ 
ple (m = 0 ,m— 2),(m = 1 ,m — 3), etc.) and are observed as doublets in the spec¬ 
trum, which is to say have almost identical dispersion curves in the case of a small 
refractive index contrast. The modes at m = 0 have an additional split due to the in¬ 
dependence of radial E p and angular components, thus triplets should be observed 
in the spectrum. 

It is convenient to distinguish modes with the dominant E p and E$ components 
and call them TM and TE-like modes, respectively. In the case of m = 0 we have the 
full analogy with slab waveguides: the TM-like (2.65b) and TE-like (2.65a) modes 
in the cylindrical case are described with similar equations as used for TM (2.30) and 
TE (2.29) modes in slab waveguides. The behaviour of TE and TM-like modes in 
cylindrical case resembles the behaviour of TE and TM modes in slab waveguides, 
which is to say TE-like group of modes are relatively unmoved to changes in the 
refractive index of external medium. Hence, approximately 50% of modes can be 
obtained correctly with the scalar model even in the case of high refractive index 
contrast. 

It is interesting to note that the so-called LP modes occur at m — 0, where 
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E p and components are decoupled. Thus, one of the components can be chosen 
independently of the other. Hence, when the refractive index contrast between the 
core and cladding is small and two modes coincide, the LP (linearly polarized) or 
CP (circularly polarized ) modes can be contracted by superposing the E p and E^ 
field components. However, this procedure is no longer correct if the refractive index 


contrast is high. The modes 



and 


E P 

0 


become separated and can no longer 


be used as a basis for construction of LP or CP modes. The name LP modes can be in 
a way misleading if the weakly guided approximation is not valid. In the general case 
at m — 0 we have a set modes, coming in pairs, with a slightly different propagation 
constants and distinct alignment of the vector E, either radial or tangential. Thus 


the energy can be coupled separately into 



or 


E P 

0 


modes. 


The split between the two LP modes can be viewed as a manifestation of light’s 
intrinsic degree of freedom. The observed split of degenerate states resembles the 
Zeeman effect where the electron energy levels are spit due to the electrons intrinsic 
degree of freedom. 

We also note that the presented dispersion curves showed that the modes posi¬ 
tioned close to the cutoff region are the most sensitive to the refractive index changes 
in the surrounding medium. The largest split between the TE-like and TM-like 
modes is also observed at some small distance away from the cutoff region, thus we 
can particularly target those modes close to the cutoff region to achieve the largest 
sensor sensitivity. 




Chapter 3 


Modeling of tilted Bragg grating 
(TFBG) structures 


The goal of this chapter is to obtain the exact solution to the problem of light 
propagating through a cylindrical waveguide with a tilted periodic structure inscribed 
along its longitudinal axis. 

3.1 Derivation 

The problem geometry is shown in Figure 3.1. First, let us consider the scalar 
Helmholtz equation (2.14): 


(V 2 + k 2 )u(r) — f(f), (3.1) 

here /(r) is the excitation term. Assuming that the excitation is a harmonic function 
f(r,t) rs j f (r)e Jujt , the solution for a linear system should also have a harmonic 
function u(f,t ) ~ u(r)e Jujt . 

In a cylindrical coordinate system the equation (3.1) can be written in the following 
form: 

(^d p {pd p ) + + d 2 z + k 2 0 n 2 ( Pl 0, z)^ u(p, 0, z) = f(p, 0, z), (3.2) 
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Perturbation in refractive index 
(Tilted grating) 



Cladding 


Coating layer 


Core 


Figure 3.1: The schematic representation of the problem geometry. 


or 

Lu(p, (j),z) — f (p, 0, z), (3.3) 

We should note that p, (j) and z coordinates are coupled through the refractive index 
n = n(p, 0, z ) function, dependent on p, 0 and z coordinates. Thus we have a coupled 
partial differential equation. For a non-tilted grating the coupling would only occur 
between z and p coordinates, but for a tilted grating all three coordinates ( z, p and 
0 ) are coupled. 

If the coupling is relatively small we can use perturbation theory. First we find 
solutions for the case of unperturbed decoupled equation, and then, using these so¬ 
lutions as a basis we solve the coupled problem in terms of this basis. 

It is in our interest to keep the coupled term as small as possible. Considering 
that the perturbation along the z axis is significantly smaller than the perturbation 
along the radial direction, the following decomposition of the refractive index profile 
can be written: 

n 2 (p,(j),z) ~n„(p) + A(p,0,F), (3.4) 

Now the coupling term A(p, 0, z ) is as small as possible. 

The next step is to find the basis functions e^r) in terms of which we can represent 
the solution of the coupled problem. There are several possible ways the basis function 
efc(r) can be chosen, however it is convenient to use the natural basis of the radial 
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part L p of the L operator (3.3). Such basis functions should satisfy the orthogonality 
condition (2.16) with respect to the weighting functions. As well, they should form 
the complete set of functions, spanning the solution space. 

Before we continue let us split the operator L, defined in (3.2), (3.3), into the 
radial and longitudinal parts. With this in mind, we search for a solution in the form 
of the series: 

u( l >,4,,z) = Y 1 c m (t>,zy m4 '- (3.5) 

m 

From now on we will use the upper indices to referencing the functions dependent on 
the angular coordinate. 

Ensuring that the function u(p,<p,z) at 0 = 0 has the same value as at 0 = 27r , 
we conclude that m = ±1, 2, 3,... for 0 e [0, 27 t]. 

Plugging the series expansion (3.5) into the Helmholtz equation (3.2) and consid¬ 
ering (3.4) we get: 


{~ p d p(p d p) -~^r + k o n2 o(p) + d l + fc o A (/h0 , z )) c?n (/h z ) ejm(p = f(p,0, z )i 


or 


Y ( L 7 + + *4Na 4>, Z)) c ™(p, z)e** = f(p, 4>, z). 


(3.6) 

(3.7) 


Here the operator acting on p is denoted as L™, and it is the radial part of the L 
operator: 

L™ = -d p (pd p )-^ + ky o (p). (3.8) 

p p 2 

Finally the basis functions can be found by solving the eigenvalue problem: 


L^(P) = Ket(p). (3.9) 

Here e™(p) are the eigenfunctions, spanning the solution space, which we are going 
to use as the basis function to solve the coupled problem (3.2). The eigenvalues are 
defined as \™. 

Let us assume that the basis functions e™(p) were successfully found, and we can 
continue solving the initial coupled equation (3.2). At this point the coupling appears 
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through c m (p,z) and A(p, 0, z) terms in equations (3.5) and (3.6). The coordinates 
z and p can be decoupled by expanding c m (p, z) term into a series over the basis 
functions e™(p), similarly to what we did in (3.5): 

UAz) = £<TU)e?(p)- (3.10) 

k 

Upon inserting the series (3.10) into (3.7) we get: 

E E ( L ™ + </>, z)) cnzKip)^ = f( P , <p, z). (3.H) 

k m 

The main advantage of our basis now becomes evident. Let us consider the eigenvalue 
equation (3.9) and replace L™e™{p ) with \™e™(p) : 

*5 A (P, </>, ^)) CT(z)eT(py m * = HP, z). (3,12) 

k in 

We are left with the system of coupled ordinary differential equations of only one 
variable z. 

Let us simplify the following derivation by introducing 0 variable: 

C(A <t>) = ejWt (3.13) 

thus the equation (3.12) becomes: 

E E ( A ? +*)) crwrtp. <« =/u </>, z). p.u) 

k in 

The upper and lower indices are used to refer to a particular mode family and a 
particular mode in the given family, respectively. 

Now let us consider the orthogonality properties of the 0) functions: 



P^i^*kdpd4 


p 2 n /* oo 

/ / p(e" i (p)e>''*)(eZ(p)e-> m *)dpd^ 

Jo Jo 

poo p2n 

I peUp)e?(p)dp I P in *e- lm *d4’ 
iTSa 2?rf mn , (3.15) 
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where the norming factor 7 ™ is defined as: 

poo 

IT =< >= / penp)eT(p)dp. (3.16) 

Jo 

Here we have considered the orthogonality relation of harmonic functions and the 
basis functions (2.107). 

The above expression (3.15) is non zero if and only if i k and m ^ n simultane¬ 
ously, otherwise it is zero, thus we have proved that 0 ™(p, 0 ) functions are orthogonal. 
Using this property let us take the scalar product of both sides of equation (3.14): 



pr^Tdpdcfr 


Cl\z) 


+ 


X J 


k.m L ' 


<T(*) = 


J Pip*if(p,<l>,z)dpd(l>, 


(3.17) 


here 0 *”(p, 0 ) = e"(p)e -jn< ^ is the complex conjugate of 0 "(p, 0 ) function. 

Finally, considering the orthogonality relation we get the system of differential 
equations: 

(A” + <®C?(z) + £ £ [V.rM] C?(z) = F?(z), (3.18) 

k m 

here 

-| -| p2ir poo 

Z7r 7 i Jo Jo 

-t -t p 2n poo 

KW = T— / / 

Ti Jo Jo 

poo 

= / pe?(p)eUp)dp. (3.19) 

Jo 

Once the system (3.18) is solved for C^\z) functions, the final solution can be con¬ 
structed with help of (3.5) and (3.10): 

u(p, <M = £ £ <T (3.20) 

m k 

The system of equations (3.18,3.19) allow us to model the process of light prop¬ 
agating through a cylindrical waveguide, with a tilted grating inscribed along its 
longitudinal axis. 







Chapter 3. Modeling of tilted Bragg grating (TFBG) structures 


57 


Let us remark here that equations in system (3.18) are enumerated by an unique 
pair of indices. In analogy with quantum mechanics we can refer to the upper index 
as the orbital quantum number, pointing to a particular mode family with a unique 
angular momentum; and refer to the lower index as the main quantum number, 
determining a particular mode in the given family of modes. 

In the most general case of vectorial modes (2.21) the eigenvalue equation is 
defined as follows: 

(dP p + f p + (hie)'d p -^ + ek* + Qne)'' -if \ (E\ (E p \ 

v j: f +?(w <+h - +<) w W 

(3.21) 

or 

KKXp) = AreT(p). (3.22) 

Here e = e(p) is considered to be independent of z and 0 variables. The dependence 
is included as a small perturbation in the A (p, 0, z) function on the next step as it is 
described above. The following derivation is identical except that now we shall con- 

( u m \ 

] instead of the scalar basis function 

v kJ 

e™(p). As a result, the form of equations (3.18) is preserved, but the matrix coeffi¬ 
cients are computed differently. The scalar product has to be modified in accordance 
with (2.105): 

POO 

7 T =< >= / (6 (pK(pK(p) + <(pK(p)) dp. (3.23) 

Jo 

The remaining steps are analogous to the scalar case. 

In this section we showed that the initial problem, defined by the partial differen¬ 
tial equation (2.14) or (2.21), with coupling along all the three coordinates z, 0 and 
p, can be successfully reduced to a system of ordinary differential equations, coupled 
only along the z axis. The coupling was introduced through the matrix elements 
[M^ m (z)], varying along the 0 axis and defined by the grating profile A (p, 0, z). The 
basis functions ef(p) are defined by the equation (2.106) and depend on the waveg¬ 
uide radial profile. We will determine the matrix elements and basis functions in the 
following sections. 
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3.2 The matrix elements 

In this section we compute matrix elements [M^(z)\, defined in equation (3.18), 
for a tilted fibre Bragg grating. Next, energy transfer from the fundamental core 
mode into a different set of cladding modes can be determined with the help of the 
coupled-mode theory (CMT). 

The approach based on the coupled-mode theory (CMT) with application to the 
tilted fibre grating was initially introduced by Erdogan and Sipe [53], following by a 
series of publications presenting numerical results [54, 9]. Alternatively, the method 
based on antenna theory and equivalence theorem [55], called the volume current 
method (VCM), can be implemented as described in [56, 57, 58], where the scattered 
field is represented by radiation of an array of elementary current dipole located 
at the index perturbation. Thus, instead of approaching the problem from a pure 
mathematical point of view, the physical analogy between a perturbation in refractive 
index and an antenna can be established. The comparison between CMT and VCM 
methods was presented in [59]. 

The tilted fibre Bragg grating is shown schematically in Figure 3.2. 



Figure 3.2: Illustration of a tilted fibre grating 

The grating is specially modulated by a sine function written with fringe planes 
that are blazed with respect to the optical axis. We will see later that the tilt of the 
grating is the key element for selective and polarization-dependent light coupling. 
Assuming a uniform and periodic grating along the z-axis, the refractive index 
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Figure 3.3: Illustration of the refractive index perturbation along the £-axis 


perturbation in the core can be represented in the following form: 

<t>, z) = n a (p ) + A n(p, 0, z), 


as shown in Figure 3.3. 


n 2 = n 2 + n 0 An + (An) 2 
~ n 2 + n a An 


according to (3.24) 


thus 


n 2 = n 2 + A 


(3.24) 


(3.25) 


(3.26) 


A (p, 4>, z ) = n 0 (p)An(p, 0, z) = 

= n 0 (p)S(p) cos(Kgz') (3.27) 

Let us assume that the z'-axis is the grating axis, tilted at the angle 6 g with 
respect to the fibre axis z. The fringe planes of the grating are written parallel to the 
y-axis, as shown in Figure 3.4. 

The coordinate system Oxyz can be introduced in such a way that the z'-axis 
transforms into the z-axis by the coordinate system rotation above the y- axis. 

Considering the rotation about the r/-axis 


z ' = xsin(6 l 9 ) — zcos(dg) 


(3.28) 
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Figure 3.4: a) A schematic representation of a tilted Bragg grating inside the fibre, 
b) projection onto the zx and c) xy planes, d) superposition of normal gratings taken 
at various </> angles. 


the equation (3.27) can be rewritten in the following form: 


A(z, p, (f> ) = n 0 (p)S(p) cos (K g xsm(9 g ) - K g z cos (9 g )) = 

= n 0 {p)8(p) cos (~K z z + K t {<j>)p) = 

= ^ e - jK * z e? KtWp + c.c. (3.29) 

Here we considered that x = p cos(0) as shown in Figure 3.4, and 

K = *L 
9 A/ 

K z /i^cos(0g), 

K t {4>) = K x cos{4>) = K g sm(9 g ) sin(^). (3.30) 


It should be noted that K t (<j>) depends on the angle 0, as shown in Figure 3.4 d), 
and reaches the maximum value along the auaxis: 
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and the minimum value along the y-axis: 

K, mn = K t (4, = 0,U = 0, (3.32) 

hence, the grating period for the transverse grating depends on the </> angle A p (</>) G 
w 2?r m oo). 

Kg Sin(0g) ’ y 

Therefore, mathematically such a grating may be described by superposition of 
infinitely many individual gratings written perpendicular to the fibre axis. 

Finally, considering (3.28) we can rewrite the matrix elements (3.18) in the fol¬ 
lowing form: 

i i f 2jr r°° 

«rw = r- / kl^P,<P,^(pK'(p)e i{m - n) *dpdd,= 

27r 7 i Jo Jo 


]_^_ e -jK z z 


471 ryn 


r-2 n 


e jKt(rp)p e j(m-n 


/ %(p)<Z(p)no(p)6(p) 

'o v^o 

i k 2 r°° 

J (p)n,MHp)<r m "(p)dp + c.c. 


dp + c.c. 

(3.33) 

Here cr mn (p) is the function dependent only on p (a cylindrically symmetric function): 


r 

a mn {p ) = e 

Jo 

r 2n 


jK t {<j>)p e j(rn-n 


= 27r(-l) fc J fc (7p), 


(3.34) 


here Jk{ip) - axe Bessel function of first kind, k = m — n, and 7 = K g sin(0 9 ). 

Hence, in the case of a tilted grating the orthogonality between the modes is 
broken due to the n 0 (p) Jk{ip) term, which arises via the perturbation caused by the 
tilted grating. 

We are mainly interested in the energy transfer from the core modes into the 
cladding modes, considering the weighted function Jk{ip) and the fact that the core 
mode is mainly confined inside the core, and the grating perturbation also exists 
inside the core, we conclude that the coupling is possible only to a limited number of 
higher order azimuthal mode families, with azimuthal number m < 12 , as shown in 
Figure 3.5. 
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core 


waveguide radius, p (p m) 


Figure 3.5: The weighted function J m ( 7 P) caused by the grating tilt, calculated 
for 9 = 0° and 9 = 10° grating tilt angles. The grating period is assumed to be 
A g = 0.6 /im, thus for 9 = 0° we get 7 = K g sin(0) = 0, and for 9 = 10° we get 
7 = 1.81 pm' 1 . The number m = m\ — m 2 is the difference between azimuthal order 
of the first m\ and the second m 2 mode families. 
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If the grating in non-tilted, the equation (3.34) is reduced to: 

a mn (p ) = / = 2v r8 mn , (3.35) 

J o 

thus, the modes belonging to different families are mutually orthogonal, and hence, 
the energy transfer between such modes is impossible. 

We conclude that the expression for coupling coefficients (3.18) is reduced to 
expression describing orthogonality between radial components of modes, except that 
we have an extra weighted function n 0 ( y p)a mn {p). The extra weighted function not 
only breaks orthogonality between modes within a particular family, but also between 
modes belonging to different families of modes with different azimuthal numbers m. 
We note that the radial components of modes belonging to different families were not 
orthogonal initially, as shown in Figure 2.16, but the complete modes with the angular 
part included, in the case of a non-tilted grating, were orthogonal due to (3.35). 

According to orthogonality relation (2.105), for the vectorial case the term 
ef (p)e^(p) in equation (3.33) should be replaced with: 

7(p>T(p) = 6 (p)p<(p)<(p) + pv?(p)v?(p), (3.36) 

are the eigenvectors of non-perturbed problem (2.21). 

Whereas in the scalar case, in accordance with (2.111), we have: 

%(p)<%{p) = pe?{p)e?(p), (3-37) 

with eigenvectors e™(p) for the scalar eigenvalue problem (2.16). 

Finally, considering the complex conjugate part c.c. , and assuming that all the 
coupling coefficients are known, we can rewrite equation (3.33) in the simple form: 

MTU) = cos (K z z)[MST], (3.38) 

here the number M^ 71 defines overlap (or coupling) between two modes: (n, i ) and 
(m, k), as it follows from equation (3.33). All such coupling coefficients are assembled 
into the matrix 


here e)"(p) = 


u 
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Let us compute coupling coefficients C™ between the core and cladding modes 
(here m is the mode family with azimuthal number m, and k is the mode order in the 
family). The results are presented in Figures 3.6, 3.7 for 2° degree, Figures 3.8, 3.9 
for 4° degree and in Figures 3.10, 3.11, 3.12 for 10° degree gratings. We note that it 
is sufficient to consider only 4, 5 and 9 families of modes for 2°, 4° and 10° degree 
tilted gratings, respectively. 

In the following section we apply the coupled mode theory to compute the spectral 
response of the gratings. The computed coupling coefficients C™ are going to used 
as input parameters. 
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Figure 3.6: The 2° degree grating assisted coupling coefficients between the core and 
cladding modes, An = 10 -4 . 
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Figure 3.7: The 2° degree grating assisted coupling coefficients between the core and 
cladding modes, An = 10 -4 . 
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Figure 3.8: The 4° degree grating assisted coupling coefficients between the core and 
cladding modes, A n = 10 -4 . 
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Figure 3.9: The 4° degree grating assisted coupling coefficients between the core and 
cladding modes, An = 10 -4 . 
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0 = 10° grating 



Figure 3.10: The 10° degree grating assisted coupling coefficients between the core 
and cladding modes, An = 10~ 4 . 
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Figure 3.11: The 10° degree grating assisted coupling coefficients between the core 
and cladding modes, An = 10~ 4 . 
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Figure 3.12: The 10° degree grating assisted coupling coefficients between the 
and cladding modes, An = 10” 4 . 
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3.3 Coupled-mode theory, 

the two modes approximation 

The coupled mode theory was employed in electromagnetics in the early 1950’s, 
and was initially applied to microwave travelling wave devices [60, 61, 62], In the 
early 1970’s the coupled mode theory had been successfully applied to the modeling 
of various guided wave optical systems. The noticeable papers were published by 
Kogelnik [63, 64], Snyder [65], Yariv [66], Marcuse [67] and others. The application 
to parallel waveguides can be found in [68, 69, 70]. 

In this section we continue our consideration of the cylindrical waveguide with 
a periodic grating inscribed in its core. The grating allows the transfer of energy, 
initially confined in the fibres core, outside the core into the cladding modes. It is 
also known that in our particular case the contra-directional phase matching condition 
occurs, as we will explain it later. 

Let us start with equation (3.18). Assuming that there is no external excitation 
we get: 

(A” + <®C?(z) + Y E KTM] c rw = 0. (3.39) 

k m 

For convenience we enumerate modes and coefficients with single index, e.g. each pair 
of indices (z, n ) we will be enumerated with a single index a, thus we can rewrite (3.18) 
as 

(A* + d 2 z )C a (z) = -J2 Wen(z)] C 7 (z). (3.40) 

7 

If the perturbation can be neglected, i.e. the right hand side of (3.40) can be 
canceled [M ai (z)\ = 0, and hence the solutions would be C a (z) = A a e zy ^* z , where 
A a are some constants. If the perturbation is switched on, the constants should be 
replaced with a slowly varying functions A a = A a (z) along the z coordinate [63, 66]. 

C a (z) = A a (z)e^ = 

= A a (z)e^ z , (3.41) 

here we denoted {3 a := called the propagation constants. 
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The final solution to (3.20) can be written in the following series form: 


u(p, 0, z) = A k( z ) e k(p ) 

m k 

Inserting (3.41) into (3.40), and considering that 


e jm<j> e ip a z 


(3.42) 


d 2 z C a (z) = d 2 z A a {z)e J ^ z = e 3 ^ z [-\ a + j2/3 a d z + d 2 z ]A a (z), 


(3.43) 


we obtain 

e j ^ z \j2M z + d 2 z \A a (z) = [M ai (z)} Ay(z)e^ z . (3.44) 

7 

The matrix coefficients M ai (z) are defined in accordance with (3.38): 

[M ai (z)\ = [M ai \ cos(K z z) (3.45) 


The equation (3.44) can be rewritten in the following form: 

\j2/3 a d z + d 2 z ]A a (z) = - ^ [Meryl A z( z ) cos (K z z)e 3 ^~ pa)z . 


(3.46) 


Now let us consider the oscillating terms: 


cos [K z z)eP^ 


J (P~f Pa)z — _ f gjK z Z _|_ g jK z Z ^ 7 Pa) 


= A 


jj (4" /tv T Kz ) _|_ (3~f pa K z ) 


(3.47) 


The matrix terms oscillate periodically and rapidly unless the phase is close to 
zero, in which case the power coupled between modes accumulates coherently and 
gives rise to a significant power exchange between the modes. Thus we can neglect all 
the matrix terms except a few terms that have a small or zero phase A0 = 0 under 
the exponential functions [71, 66, 50]. 

Considering (3.47) and assuming that the energy couples from the core mode 
with the propagating constant f3 a , to one of the cladding modes with 0 7 propagating 
constant, there are two possible phase matching conditions: 


A cj) = K z + 0 7 - /3 a = 0, 
A0 = K z + fd a - Ay = 0, 


(3.48) 
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called co-directional and contra-directional coupling, respectively. 

From the experiments described in the following sections we know that in onr case 
only the contra-directional coupling occurs, thus we limit ourselves to this particular 
case. Figure 3.13 shows the phase matching condition, or momentum diagram, of 
contra-directional coupling. 


◄- 

-► 

-W- fa fa 

s 

Figure 3.13: The momentum diagram of contra-directional coupling (here fa corre¬ 
sponds to the core mode and fa to the cladding mode). 

Now let us make a few extra assumptions. First we assume that instead of consid¬ 
ering coupling between all the modes at once, we can consider only two modes, the 
core mode and one of the cladding modes. This assumption is useful to understand the 
basic properties of energy transfer from the core mode to a particular cladding mode. 
Once the energy is transfered to a particular cladding mode it can not be recoupled 
to an other cladding modes due to the significant phase mismatch. We also assume 
that the amplitudes A a (z ) in (3.41) are “slowly” varying, i.e. d 2 A a (z ) -C f3 a d z A a (z), 
this is the so-called slowly varying amplitude approximation [63, 66]. The resulting 
equations (3.46) will takes the following form: 

j2fad z A l (z) = -Ke~ jSz A 2 (z), 

j2fad z A 2 (z) = - K *e iSx A 1 (z). (3.49) 

Here: 

6 = —K z + fa — fa —> 0 is the phase mismatch. We should keep in mind that if the 
propagation constant fa is positive then, due to the contra-directional coupling the 
fa constant should be taken with the minus sign. 

« = [Mi 2 ] = [M 21 ]* =< > is the coupling constant, and can be computed in 

accordance with equation (3.18). 
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We note, by looking at the equation (3.47), that the off-diagonal matrix elements 
have phases approaching zero in the case when the phase matching condition is sat¬ 
isfied, thus indeed we have condition that allows for coherent power accumulation, 
which leads to the power redistribution between the modes. 

Let us solve equations (3.49) assuming that the phase matching condition is satis¬ 
fied 5 = 0, and the proper boundary conditions are given. Considering that the core 
mode with the amplitude Ai(z) to be incident at z = 0 on the perturbation region 
z G [0, L] we can set 7^(0) = 1. The cladding mode A 2 (z) is “generated” by the per¬ 
turbation, hence A 2 (z) should be set to zero at the end of the perturbation region 
A 2 (L ) = 0. The solution of (3.49) (for the case of 5 = 0) is given by equations (3.50) 
and mode power of the incident and scattered waves are shown in Figure 3.14. 



(3.50) 



ML) 


n 


Mi(0)| 2 



Figure 3.14: The transfer of power between the incident core mode A\(z) and back- 
scatters cladding mode A 2 (z) in the case of contra-directional coupling. Here ^ = 2.4 


Now let us consider the amplitudes A\ and A 2 at the perturbation boundary. 
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Setting as previously 7Li(0) = 1 and A 2 (L) = 0 we get for the general case of non zero 
phase detuning (5^0) the following solution: 


Here 


ML) 

H 2 (0) 


a cosh(L|) — i<5sinh(L|) 
-f sinh(Lf) 
crcosh(L|) — iS sinh(L|)' 


(3.51) 



^ — —kz + A — A; 


(3.52) 


here L is the perturbation length, A and A are the propagation constants of the core 
and cladding modes, respectively, as shown in Figure 3.13. 

We have a particular interest in the energy loss due to the coupling from the 
core to the cladding modes, as this energy loss can be observed experimentally by 
measuring transmission spectra of the optical system: 


rj(x) 


\ML,x) l 2 

Ai(o, x)\ 2 

_7 2 (A_ 

7 2 (x) cosh 2 (C^^) + x 2 sinh 2 (C' 2 ^) ’ 


(3.53) 


here C is the coupling parameter defined as C := La, x | is the phase detuning 
parameter, £ := k/\/AA and 'y(x) := | = a/ 1 — x 2 . 

The power loss r](x) function of the core mode due to the coupling to the cladding 
mode is shown in Figure 3.15 for various coupling parameters C = La and the phase 
mismatch x — 

In the case of coupling to the lossy mode we note that the peak becomes broader 
and the base line is shifted, as shown in Figure 3.16. 
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Figure 3.15: The power loss of the core mode as a function of the phase mismatch, 
for various coupling parameters C = La. 




wavelength, A (pm) 

b) 


Figure 3.16: The power loss of the core mode as a function of the phase mismatch, for 
the fibre waveguide coated with a thin him (h = 100 nm) made of lossy material (blue 
curves) a) N eff = 1.3661 - j6.1165 x 10~ 4 , and b) N eff = 1.3669-^2.2791 x 10" 5 . 
The red curves, used as a reference, corresponds to the non-lossy material ( 
N eff = 1.3661 ). 
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3.4 Coupling between the core mode and many 
cladding modes in the TFBG 

In this section we apply the coupled mode theory to our problem of interest: the 
tilted fibre Bragg grating (TFBG) inscribed inside the core of a standard telecommu¬ 
nication fibre SMF-28. 

The phase matching condition in such a case is schematically shown in Figure 3.17, 
where the core mode is coupled to a number of cladding modes. 



Figure 3.17: The phase matching condition in TFBG 

In the case of TFBG the energy confined inside the fibres core is coupled from 
the core mode into a multiplicity of cladding modes, as shown schematically in Fig¬ 
ure 3.18. The grating-assisted coupling, for the problem of interest, is only possible 
for the core and cladding modes propagating in opposite directions. It should also be 
noted that for the given geometry and operational wavelength the energy coupling 
between different cladding modes is prohibited due to the significant phase mismatch. 

If a spectrum of TFBG grating was obtained from an experiment, the grating 
wavenumber Kq can be expressed in terms of the measured value of the Bragg res- 
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Figure 3.18: The schematic ilustration of energy transfer from the core mode into the 
cladding modes in the TFBG. The various potential barriers correspond to different 
azimuthal numbers m. 


onance position A b, by assuming that the forward propagating wave in the core is 
coupled to the backward propagating wave in the core: 

97r 

K g = 2 fa = 2 N b — . (3.54) 

The phase detuning condition condition (3.48) takes the following form: 

Akj(X) = 


with 

Pb(X) 

(A) 

Here /3b and [3j are propagation wave numbers of the core mode, with the effective 
refractive index Nb, and the j'-th cladding mode with the corresponding effective 


K G -{p B { A)+ j 9 i (A)) = 
N b Ns + NjiX) 


2tt 2 


A 


B 


A 


(3.55) 


2n 

= N B k Q = N b —, 
= N j (X)k 0 = N j (X) 


2tt 


(3.56) 
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refractive index Nj, computed at the operational wavelength A. We have also assumed 
that Np( A) ~ constant approximately does not depend on the A, which is true if the 
radius of the core is significantly smaller then the radius of the cladding. 

The Bragg Condition is obtained by assuming that there is no phase detuning 
between the modes, i.e. A kj ~ 0. Hence the position of the j-tli resonance is defined 
by the following equation: 


A 


resj 


1 

2 



N b ) 


A B 


(3.57) 


Assuming the operational wavelength A and considering (3.55) we can determine 
the phase detuning for the j'-th resonance in the proximity to Ay 


5j(Xj — A) = 2tt 2 -A- 


N b Ns + Nji A) 


A 


B 


A,- 


-2tt 2 


N b N b + Nj(X) 


A 


B 


A 


— 2tt(N b + Nj(X)) ( —— + - 
2 tt(Nb + N,( A)) ( 


(3.58) 



wavelength, A (p m) 


Figure 3.19: The resonances of different families of modes, computed in accordance 
with (3.57). The peaks arising due to the coupling between the core mode and m — 0 
family of modes, neglecting the coupling to higher family of modes. Here the grating 
length L = 10 mm and the coupling constant C — 2 • 10 -4 . 
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Using the two mode approximation and considering the dephasing expression (3.58) 
the resonances of a particular family of modes, let us say for m = 0, can be plotted 
along the wavelength of operation A, as shown in Figure 3.19. For a single family of 
modes the resonance peaks have almost no overlap (Figure 3.19), thus a resonance 
can be approximately computed with the two mode approximation independent of 
the other resonances. In other words, the energy coupling between the core mode 
and a particular cladding mode can be computed without considering other cladding 
resonances of the same family of modes. 

However, in the case of TFBG, many resonances, corresponding to various families 
of modes, are present in close proximity to each other. In Figure 3.19, resonances of 
the first 11 families of modes are plotted along with the peaks computed for m = 0 
family. Although the two mode analysis is useful to study basic properties of the 
problem, it can not be applied to our problem of interest due to the overlap between 
the resonances. 

Going back to the initial coupled modes equation (3.46) and, as previously, ne¬ 
glecting the rapidly oscillating terms, not contributing to the change of amplitudes, 
we obtain the following system of coupled equations: 

N r 

d z A core {z) = 

k =1 ^ 

d z Ai‘ d (z) = (3.59) 

Here Ck =< ^ core \A\^ lad > is the coupling constant between the core and the k-th 
cladding mode, as was discussed in the previous section; N is the number of coupled 
modes. The values of Sk( A), /3% lad (\), ^(A) and Ck( A) have to be computed at a 
particular operational wavelength A. 

The analytical solution, similar to (3.50), can no longer be applied to the system 
of equations (3.59). Instead this system of coupled nonlinear differential equations 
can be solved numerically. 

In the case of co-directional coupling the initial value problem (IVP) is considered 
A core (0) = 1 and A^ ad ( 0) = 0, hence the system of equations can be directly inte¬ 
grated. Unfortunately, in our case of contra-directional coupling the boundary value 
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problem (BVP) A core ( 0) = 1 and A% ad (L) = 0 has to be considered (as the bound¬ 
ary conditions are known at the opposite sides of the interval), and this significantly 
complicates the numerical routine. Likely the problem (3.59) is linear with respect to 
the ^ variable, hence we can start propagation from the opposite end z = L assuming 
that A c ^ ad {L) = 0 and setting amplitude of the core mode at z = L to some arbitrary 
constant C to be determined later: A core (L) = C. Next the system of equations can 
be integrated as an initial value problem. Once the solution v4 core (0) at z = 0 is known 
we can renormalise the solution to ensure that A core (0 ) = 1. The procedure is shown 
in Figure 3.20. 
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Figure 3.20: Transforming the boundary value problem into the initial value prob¬ 
lem. a) solving IVP by assuming A core (L) = 1 and A clad (L ) = 0, b) renormalized the 
solution by setting A core ( 0) = 1. 


The initial value problem is next solved with the help of the standard fourth 
order Runge-Kutta method. The result for the core mode coupled to the two closely 
positioned cladding mode resonances is shown in Figure 3.21. 

In Figure 3.21 a) the difference between the two mode approximation, when the 
coupling to cladding modes is considered independently (the green and the blue 
curves), and the exact approach (the red curve) is clearly seen. It is also interesting 
to notice that the energy from the cladding mode Figure 3.21 c) (the blue curve) can 
be recoupled back to the core mode (the red curve). 
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Transmission Loss (the core mode) Reflected Energy (the cladding modes) 




Figure 3.21: Coupling between the core mode and two cladding modes. The grat¬ 
ing length L — 10 mm and the coupling constant are C\ — 2 • 10 4 and C 2 — 1 ■ 10 4 
between the core and the two cladding modes. Here Ai(z) is the amplitude of the 
forward propagating core mode (the red line), A 2 (z ) and A 3 (z ) are amplitude of the 
backward propagating cladding modes (the blue and the green lines). 


Finally, the coupling problem for the TFBG structure of our interest can be solved 
as follows: 

1. We introduce a grid vector Aj spanning the operation range of the sensor, i.e. 
Aj G [A m i n , A m ax], where at each point Aj the system of coupled equations (3.59) 
has to be solved. We note that a typical peak width is about 0.2 nm, where 
the operational range is about 100 nm. Thus, at least 1000 points Aj have to 
be considered to observe the resonances. Considering the necessity of solving 
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the system of coupled differential equations (3.59) at each point, the described 
approach can be extremely time consuming. We can overcome this difficulty by 
introducing a nonuniform grid. The positions of resonances is known (3.57): 

Ke,, = \ (l + A b, (3.60) 

hence we can introduce a finer mesh in the vicinity of resonances, as shown in 
Figure 3.22. 


non-uniform grid 
resonances 



wavelength, (p m) 


Figure 3.22: The non-uniform grid with a finer mesh in the vicinity of resonances. 


2. Next the detuning parameter Sj for each resonance can be computed in accor¬ 
dance with (3.58): 

Sj( A) ~ 2 7r(N B + Nj(\)) , (3.61) 

3. Unfortunately the dispersion of modes can not be neglected, as can be seen 
from Figure 3.23. Hence, we have to find modes and coupling constants at each 
point A j of the grid, or at least we can split the interval [A m j n , A max ] into a set of 
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smaller subintervals and compute resonances and coupling coefficients for each 
of them, considering them to be constant within the subintervals. 


X = 1.5 ((.i m) 
A = 1.6 (p m) 



wavelength, (n. m) 


Figure 3.23: The dispersion of modes. The resonances of m = 0 family of modes are 
computed at Ai = 1.5 /im and A 2 = 1.6 /im operational wavelengths. 


4. Finally the system of coupled nonlinear differential equations (3.59) can be 
solved independently at each point A j of the grid: 


d z A core (z, Xj) = E * oSla ' Aj), 

k=l 2 ' j ^ 

dX^Ai) = (3.62) 

^Pk \Pj) 


The number of coupled modes N can be limited to 15 modes by choosing only 
the modes with the smallest phase mismatch Sk(Xj)- The remaining modes can 
be neglected due to a weak coupling. For each A j the choice of N coupled modes 
should be reconsidered. 


The theoretically computed transmission spectra along with the experimentally 
measured spectra are shown in Figures 3.24 and 3.25 for the 2° degree grating, in 
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Figures 3.27 and 3.28 for the 4° degree grating and in Figures 3.31 and 3.32 for the 
10° degree grating. A more detailed information can be obtained by zooming in to 
particular resonances, which are shown in Figures 3.26 for the 2° degree grating, in 
Figures 3.29, 3.30 for the 4° degree grating and in Figures 3.33, 3.34, 3.35 for the 10° 
degree grating. The coupling coefficients are plotted on the same figures. 

We note that the resonances have a similar structure along the whole operational 
range of the sensor. We also note the an interesting effect of peaks alternation. There 
are only two group of peaks: the peaks consisting of modes with odd azimuthal 
symmetry m — 1,3, 5..., and modes with even azimuthal symmetry m — 0, 2, 6... . In 
the spectral area where peaks are separate, these two groups of peaks are alternating. 
This is an important effect explaining some of the experimental data, to be discussed 
later. 
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: 3.24: The experimentally measured spectra of the 2° degree 1 cm long grating. 



Figure 3.25: The theoretically computed transmission spectra of 2° degree grating. 
(L = 1 cm, An = 10~ 4 ) 
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Wavelength, X (u m) 


Figure 3.26: The fine structure of the particular resonances and corresponding cou¬ 
pling coefficients of 2° degree grating. 
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: 3.27: The experimentally measured spectra of the 4° degree 1 cm long grating. 



Figure 3.28: The theoretically computed transmission spectra of 4° degree grating. 
(L = 1 cm, An = 10~ 4 ) 
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Figure 3.29: The fine structure of the particular resonances and corresponding cou¬ 
pling coefficients of 4° degree grating. 



Wavelength, X (li m) 

Figure 3.30: The fine structure of the particular resonances and corresponding cou¬ 
pling coefficients of 4° degree grating. 
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Figure 3.31: The experimentally measured spectra of the 10° degree 1 cm long grating. 



Figure 3.32: The theoretically computed transmission spectra of 10° degree grating. 
(L — 1 cm, An = 10~ 4 ) 
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Figure 3.33: The fine structure of the particular resonances and corresponding cou¬ 
pling coefficients of 10° degree grating. 



Wavelength, X (u m) 

Figure 3.34: The fine structure of the particular resonances and corresponding cou¬ 
pling coefficients of 10° degree grating. 
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Figure 3.35: The fine structure of the particular resonances and corresponding cou¬ 
pling coefficients of 10° degree grating. 

As can be seen from the presented figures the developed theoretical model provides 
results in good agreement with the experimental measurements. However, the small 
difference is present in the area of the so-called ghost modes [11] for 4° degree tilted 
grating. This area is highly populated with resonances, thus our approximation where 
only 15 coupled modes are considered might not be sufficient. It also should be noted 
that we presume an ideal grating. In reality a physical grating might be subjected 
to various unaccounted-for effects, such as a non-uniform and one sided illumination 
by the UV beam in the process of grating inscription. Nevertheless, from the point 
of view of practical application, we are interested in the spectral area with distinct 
resonances. This region is simulated in good agreement with the measurements. 
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3.5 Polarization-dependent coupling 

In this section we discuss polarization-dependent properties of the TFBG. Con¬ 
sidering that the core mode has rn — 1 azimuthal symmetry it can be orientated 
differently in the plane transverse to the propagation axis. The tilt of the grating 
planes breaks the cylindrical symmetry of the fibre, and defines the reference frame 
x — y in which it is convenient to analyze the system. Considering the mutual ori¬ 
entation of the grating planes and the incident core mode, we can expect that the 
coupling coefficients might depend on the relative orientation between the polariza¬ 
tion (or orientation) of the incident core mode and relative orientation of the grating 
planes, as shown in Figure 3.36. 



Figure 3.36: The schematic representation of the TFBG grating and the incident 
linearly polarized core mode (the E p component is shown). The core mode is rotated 
about the optical device axis by some angle r. The grating is tilted by 6 angle about 
the x axis. 

The coupling coefficients C™ 1 in accordance with (3.33) are proportional to 

Rnp)K(p)n 0 (p)6(p)a mn (p, r)dp + c.c., 


CT(r) 




o 


(3.63) 
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here r is the angle at which linearly polarized light is incident on the TFBG structure 
and <7 mn (p, t ) is the weighted function, shown in Figure 3.37: 


-2 TT 


<J mn (p,T) = 


ojK g sin(0g) sin (tp+r)p j (m-ri 


(3.64) 


9![a mn (p, x)], for m - n = 0 
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Figure 3.37: The weighted function a mn (p , r) in the fibre core, for the 4° degree tilted 
grating. 

As it shown in Figure 3.37 the weighted function a mn (p,r) corresponding to cou¬ 
pling between the modes of the same family (in our case the core mode has the 
azimuthal number m — 1 is coupled to the m = 1 cladding modes) is not affected by 
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the incident angle r of linearly polarized light, whereas coupling to different families 
of modes, with the m — n = 1,2, 3, reveals a significant angle dependence. 

Hence, considering the weighted function a mn (p 1 r), the polarization-dependent 
coupling coefficients can be computed. For instance for the 4° degree tilted grating 
the coupling coefficients are shown in Figure 3.38 for t = 0°,45°,90° polarization 
angles of the incident light. 

Now let us study the polarization dependence of particular resonances in more 
detail. As we mentioned in the previous section there are only two alternating groups 
of resonances with either odd or even azimuthal symmetry, thus it is sufficient to 
analyze two closely positioned resonances. The results for 4° and 10° degrees gratings 
are shown in Figure 3.39 and Figure 3.40, respectively. 

It is clearly seen that the resonances have a complicated inner structure. Sev¬ 
eral coupling coefficients of various polarization-dependent behaviors are bounded 
together by the system of coupled mode equations. We note that by changing the 
polarization of incident light the dominant coupling coefficient can be changed, thus 
the energy can be predominately coupled to a particular mode with a specific field 
distribution at the sensor surface. 
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icjxicf 


0 = 4° TFBG, at 0° polarization 



Figure 3.38: The coupling coefficients Cfc(r) for the 4° degree TFBG computed for 
r = 0°,45°,90° polarization angles of the incident light. 



































































































Chapter 3. Modeling of tilted Bragg grating (TFBG) structures 


98 


|C k | x 10 ' 4 4° TFBG at P = 8° polarization 





Figure 3.39: The polarization dependence of coupling coefficients corresponding to 
two resonances of the 4° degree TFBG. 
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|C | x io' 5 10° TFBG at P = 8° polarization 
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Figure 3.40: The polarization dependence of coupling coefficients corresponding to 
two resonances of the 10° degree TFBG. 
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Let us now assemble the density plot by computing the TFBG spectra at each 
angle of linear polarized light (the angle of rotation about the optical axis). First we 
compute the coupling coefficients for cladding modes of m — 0,1,2, 3,4, 5 azimuthal 
symmetry (as was previously shown the six first families of cladding modes are suffi¬ 
cient to accurately compute the spectrum of 4° degree TFBG). 



Figure 3.41: The coupling coefficients between the core and cladding modes of az¬ 
imuthal order m — 0,1, 2, 3,4, 5 computed at various angles of linearly polarized light 
incident at the 4° degree TFBG. 


The coupling coefficients of various m azimuthal numbers plotted on the same 
figure are shown in Figure 3.42. 

Applying the developed in the previous chapter method we compute the corre¬ 
sponding transmission spectra for each angle of linearly polarized light incident at 
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1.48 1.5 1.52 1.54 1.56 1.58 

Wavelength, X (jx m) 


Figure 3.42: The coupling coefficients of the 4° degree TFBG computed at various 
angles of linearly polarized light. 

the TFBG structure. The result of the computation is shown in Figure 3.43, the 
experimental measured spectra at various polarization angles are shown Figure 4.1. 


I,(dB) 
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Figure 3.43: The transmission spectra of the 4° degree TFBG, computed at various 
angles of linearly polarized light. 
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3.6 The electric field distribution at the fibre 
boundary 

In the following Chapters we discuss the possible ways of improving the TFBG 
sensor sensitivity by depositing nanoparticles on the sensor surface. The light par¬ 
ticle interaction is dependent, among other parameters, on the incident electric held 
orientation. In this section we discuss the electric held orientation at the hbre surface. 

The developed vectorial mode solver allows us to determine E p and held com¬ 
ponents for each of the modes. In Figure 3.44 and Figure 3.45 the held components 
E p and E^ are plotted along the effective refractive index of the modes with different 
azimuthal symmetries m = 0,1, 2, 3,4, 5. We note that by default the held of each 
mode is normalized to unity. 

As can be seen from Figure 3.44 and Figure 3.45 the modes with a high effective 
refractive index are almost completely bounded inside the waveguide, with almost 
zero held at the hbre surface, whereas helds with low effective refractive index modes 
leak outside the waveguide. This effect can be seen in more detail if the held is 
computed at some distance from the hbre. In Figure 3.46 the held components E p 
and E$ are computed at the distance of 1 rirri away from the hbre boundary. We 
note that these modes with low effective refractive index have long exponential tails 
penetrating deep into the surrounding medium, these are the so-called evanescent 
waves. The long exponential tail of the modes is the reason for high sensitivity of the 
sensor in the region close to the cutoff. 

The Figures 3.44, 3.45 show components of the electric held normalized to unity. 
However, the held should depend on the energy coupled to a particular mode, as each 
mode is excited with different strength. The energy (or light) coupled to a particular 
k-th mode is proportional to the corresponding coupling coefficient Ck- Thus the 
resulting held at the sensor surface can be obtained by multiplying the electric held 
components E p and E^ of the k-th mode by the Ck coupling coefficient. The result is 
shown in Figure 3.47 for various states of linearly polarized light incident on the 4° 
degree TFBG structure. 
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Figure 3.44: The field components E p and E^ computed at the fibre boundary for 
modes with m — 0,1, 2 azimuthal symmetry. 











Chapter 3. Modeling of tilted Bragg grating (TFBG) structures 


104 





Figure 3.45: The field components E p and E$ at the fibre boundary for modes 
with m = 3, 4, 5 azimuthal symmetry. 
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Figure 3.46: The field components E p and E^ computed at 1 nm distance away from 
the fibre boundary. 
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Figure 3.47: The intensity of electric field components E p and at the fibre 
boundary for 4° degree TFBG, computed at various states of linearly polarized light 
P = 0°, 45°, 90°. 
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Let us zoom in on a few particular resonances. As mentioned previously each 
resonance consists of modes of either odd or even azimuthal symmetry, as shown in 
Figure 3.48. 



Wavelength, X (g. m) 

Figure 3.48: The structure of the particular resonances of 4° degree TFBG. 

Multiplying the electric held components E p and E$ at the fibre surface by the 
corresponding coupling coefficients we get the value of electric held present at the 
hbre surface, shown in Figure 3.49. 

We conclude by stating that the electric held at the hbre surface is predominantly 
radially polarized (with electric held component normal to the hbre surface) in the 
wide spectral region, as can be inferred from Figure 3.47. However, the composition 
of resonances is complex (Figure 3.49) and consists of many modes with radial and 
tangential dominant polarization. Moreover, the composition depends on polarization 
of the incident light, as shown in Figure 3.47. Thus, by changing the core mode 
polarization, the coupling coefficient can be changed, and hence the energy couples to 
a particular mode with a specihc held distribution at the sensor surface can be changed 
as well. The spectral response of the TFBG sensor is polarization-dependent, as shown 
in Figure 3.43. The electric held at the sensor surface is polarized, with orientation 
of E dependent on the core mode polarization. Hence the spectral response of the 
TFBG sensor can be further enhanced by coating its surface with a layer, the optical 
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Wavelength, X (fi m) 



Wavelength, X (u m) 


Figure 3.49: The electric field at the fibre surface, corresponding to particular reso¬ 
nances. The linear polorizes light incident at P = 45° and P = 90° angles at the 4° 
TFBG. 


properties of which are polarization-dependent. 

In the following Chapter 4 we discuss experimental measurement techniques taking 
advantage of the polarization-dependent response of the TFBG sensor. The observed 
polarization-dependence of the TFBG spectrum and electric field at the surface of the 
sensor will lead us to the sensitivity enhancement technique proposed in the following 
chapters. 



























Chapter 4 


Experimental polarization-based 
optical sensing with application to 
TFBG sensors 


In this Chapter we proposing a polarization-based sensing method developed for 
TFBG sensors. 

Polarization-based sensing is crucial for various types of optical sensors, in partic¬ 
ular for stress analysis, plasmon-mediated sensing, and sensing of anisotropic media 
or other forms of perturbations [72, 73, 74, 75]. In the particular case of waveguide- 
type sensors (including optical fibres), it is possible to investigate the devices optical 
properties with polarized light but in general this requires very careful alignment and 
control of the input polarization. 

As we described in the previous sections, TFBG sensors reveal strong polarization- 
dependent properties [76, 3] due to the tilt of the grating planes which breaks the 
cylindrical symmetry of the fibre and strongly impacts the magnitude of the coupling 
coefficients between the incident core mode and the cladding modes excited by the 
grating [9], as is schematically shown in Figure 3.36. 

The grating tilt causes asymmetry in coupling between the cladding modes and 
the core mode, with respect to the grating tilt. In particular, for a grating tilted 
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by an angle 9 in the reference y — z plane and a linearly polarized input core mode, 
whose polarization is rotated by an arbitrary angle (j) with the reference x axis, the 
coupling to individual cladding inodes will depend strongly on both 6 and (j). 


1(90°) = | y 



1551 1552 1553 1554 1555 

Wavelength, (nm) 


Figure 4.1: A typical TFBG transmission spectrum for linearly polarized light, and 
series of spectra obtained by rotating a linear polarizer about the optical axis are 
shown as the density plot. 

The asymmetry in sensing results from the different nature of interaction between 
the cladding TM-like and TE-like modes, with the predominant radial and tangential 
orientation of the electric held at the sensor surface, from one side, and the environ¬ 
ment under test from another side. The asymmetry can be significantly enhanced 
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if the sensor surface is coated with a thin metallic hlm[l] or metal nano rods [3], 
interacting differently with the radial and tangential oriented electric fields. 

The polarization effects in optical fibres have been traditionally quantified with 
a polarization-dependent loss (PDL) parameter, provided by optical vector analyzer 
(OVA) instruments. Alternatively, the polarization effects can be studied with regards 
to linearly polarized light aligned with a specific axis of the device under test [1]. A 
typical series of spectra measured at various angles of linearly polarized light is shown 
if Figure 4.1. 

The effect of the input mode polarization can be observed by measuring optical 
transmission spectra with a polarizer inserted between the light source and the grat¬ 
ing [1], A typical series of spectra measured at various angles of linearly polarized 
incident light is shown in Figure 4.1 for a 1 cm long TFBG immersed in water. 

In the case of TFBG sensor modes with a different orientation of electric held at 
the sensor surface can be excited by rotating linear polarized light about the sensor 
axis [9]. To rotate the linearly polarized light an external polarizer is usually used, 
sometimes in combination with OVA instrument. 

In the presented work we compare the PDL based approach, which is based on 
measurements of transmission spectra along the orthogonal principal axes, with the 
measurements based on extracting linear states of polarization along a predefined axis 
either from the Jones matrix or the Stokes vector. 


4.1 The Optical Setup 

The measuring optical system was designed to acquire transmission spectra of 
tilted fibre Bragg sensor at various polarization states. 

Two alternative approaches were used. The first straightforward approach as 
shown in Figure 4.2, is based on introduction of a linear polarizer in the optical 
path [1], thus allowing the collection of spectra at various states of linear polarization. 
The second approach was based on measuring the Stokes parameters or alternatively 
the Jones matrix elements, and will be discussed later. 

The experimental setup is shown schematically in Figure 4.2 consists of a SI720 
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Micron si720 


Polarizer Controller 
PR2000 


Plating Bath 



Figure 4.2: The optical setup based on SI720 spectrophotometer and PR2000 polar¬ 
ization controller. 


spectrophotometer (Micron Optics) and a polarization controller PR2000 (JDS Uniphase). 
The SI720 spectrophotometer is based on the fast sweeping tunable laser and inte¬ 
grated photo detector. The polarization controller was set to continuously scan all 
linearly polarized light states in 80 seconds. 

The spectra over the full operational range (from 1520 to 1570 nm ) were taken 
continuously at a rate of one spectrum each 0.2 seconds, thus the spectra were acquired 
with the resolution of less than 1° degree of polarization angle. Such high resolution 
allowed us to collect a huge amount of data and, as a result, to precisely describe the 
spectral and polarization response of the TFBG-SPR sensor. 

An example of acquired spectra is shown in Figure 4.1. The spectra are represented 
in the form of a density plot. The spectra are stored in the matrix, with each row 
corresponding to a different angle of linearly polarized light. Two orthogonal states 
of polarization I x and I y were extracted at the data processing stage, and correspond 
to the measurements along axes which the maximum and minimum transmission 
spectra. 

Considering that light transmitted through a device under test is experiencing the 
insertion loss dependent on the state of incident light polarization, as a function of 
optical frequency, two alternative approaches to optical measurements are possible. 

The optical system can be fully characterized either by the Stokes parameters or the 
Jones matrix. 

The Stokes parameters and Jones matrix were acquired with help of a JDS Uniphase 
SWS-OMNI-2 system and from an optical vector analyzer (OVA) 5000 from Luna 
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Figure 4.3: The Stokes vector representation. 

Technologies, respectively. Both approaches are based on the same operational prin¬ 
ciple, except that the optical vector analyzer from Luna Technologies is also capable 
of measuring the phase delay in an optical device under test. 

An optical vector analyzer usually consists of an optical source, a polarization 
controller capable of producing one of four known polarization states and a power 
meter measuring an insertion loss. Only four measurements to determine Stokes 
parameters Sq, Si, S2, S3 are required. The Stokes parameters can be visualized with 
help of a Poincare sphere sphere as shown in Figure 4.3. First, an optical signal is 
polarized to produce one of four known polarization states a, b, c and d, described by 
the Stokes vectors So, a , *So,i» So,c and S04, and then transmitted through the device 
under test. The corresponding transmitted powers T 0ja , T ob , T 0 c and T 0jrf of the output 
signal are measured with the power meter [77, 78]. 

The transmissivity of the optical component under test can be represented by a 
4x4 Mueller matrix [M\. 

S out = [M]S in , (4.1) 

where S in = (So, Si, S2, S3) and S out = (T 0 , 7\, T 2 , T 3 ) are the input and the output 
states of polarization, respectively, represented by Stokes vectors. T 0 is the intensity 
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of transmitted light, measured with the power meter. Hence, for each state of polar¬ 
ization a, b, c, d the four transmitted intensities T 0)O , T o b , T 0>c , T 0]rf can be measured, 
and four equations can be written by multiplying the first row of the Mueller matrix 
by the input Stokes vector: 


To,a — r ^OO*S'o,a + m 01*S'l,a + ^02^2,0 + ^0 k^ 3 ,a 

To,b = ■^00*S'o,b + m 01*S'l,6 + m 02S2,fe + ^0fc*S'3,fc 

Tq,c ~ m OoSo,c + m 0lSl,c + mo2S2,c + nT' 0 kS 3 ,c 

T 0y d = '^oo*S'o,d + + m 0 2S2 )( i +'^ofc*S , 3,d 


(4.2) 


Solving the above system of equations the four matrix elements m 00 , m 0 i, m 0 2 and 
m 03 can be found, in terms of which other parameters, such as polarization-dependent 
loss (PDL), can be expressed [77]. 

Unfortunately the described technique does not allow us to gather information 
about a phase delay in an optical system. The complete information about an optical 
device can only be obtained if the instrument used can measure phase as a function 
of frequency and polarization, thus providing the complete response of a system. One 
of the possible realizations of such a device is schematically shown in Figure 4.4 [78]. 

The principle of operation is based on the so-called swept-homodyne interferome¬ 
try [78] performed for various states of polarization and fully characterizes an optical 
system by measuring the attenuation and phase delay as functions of optical fre¬ 
quency. As shown in Figure 4.4 the device consists of an unpolarized light source, 
supplying unpolarized coherent light whose optical frequency is swept continuously 
as a function of time, which is first passed through the interferometer then through a 
three-way polarization splitter coupled to individual detectors. The detection polar¬ 
ization splitter is designed to measure the light intensities Pi, P 2 , P 3 by projecting the 
light polarization state onto the three linearly independent preselected Stokes axes of 
polarization Si,S 2 ,S 3 , shown in Figure 4.3, as a function of optical frequency. The 
fourth measured P 4 is a total polarization-independent power. The optical delay is 
measured with a Michelson interferometer, where a device under test is coupled in 
one arm of the interferometer and the other arm of the interferometer is used for the 
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Figure 4.4: The principle of operation of an optical vector analyzer . 


reference. As a result a 2 x 2 Jones matrix [J] is computed at each optical frequency: 


Eout — [J]Ei : 


(4.3) 


The four elements of the Jones matrix are in general complex numbers, encoding the 
attenuation and phase delay of an optical system. 


4.2 The data processing technique 

The optical setup based on OVA implementation is shown in Figure 4.5. In spite 
of the apparent simplicity, the data analysis is significantly more complicated than in 
the case where the state of polarization is rotated mechanically. 

In this section we describe the data processing technique which allowed us to 
extract polarization-based parameters, characterizing the sensor response, from the 
data provided by an optical vector analyzer. 
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Figure 4.5: The Experimental Setup. 

4.2.1 Measurements along principal axes of an optical 
system 

In this section we review polarization-dependent loss (PDL) technique and provide 
an approach to study the system transmission spectrum along its principal axes. 

A linear optical system can be represented in terms of the 2x2 Jones matrix [J] 
connecting an incident and transmuted electric field vectors [79]: 

Eout = [J]Ei n ( 4 - 4 ) 

The transmission spectrum of a device under test usually is characterized with 
two major parameters. The polarization independent parameter called the insertion 
loss /(A) and the polarization-dependent loss (PDL) parameter, both measured as a 
function of optical frequency A. To extract these parameters from the Jones matrix, 
a hcrmitian matrix [ H ] has to be constructed first. 

[H] = PPPl (4-5) 

Next, performing eigen decomposition of [H] matrix 

[H] = [U][A\[U] T (4.6) 

we find diagonal matrix [A] containing eigenvalues p\ and P 2 of [H] matrix. 

Alternative singular value decomposition (SVD) of the Jones matrix can be com¬ 
puted [80]: 

M = IGPlTf 


(4.7) 
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with the diagonal matrix [£] containing two singular values <j\ and 02 such that 
Pi = erf and p 2 = erf, due to the fact that [H] = [J]t[J], 

The matrix [U] = (ui,u 2 ) contains two orthogonal eigenvectors (due to the prop¬ 
erties of SVD decomposition), corresponding to p\ and p 2 eigenvalues, align with 
principal axes of the system. The principal axes U\, u 2 of TFBG sensor are showed 
schematically in Figure 4.6. 



Figure 4.6: TFBG with physical x, y axes and principal system axes u\ and u 2 
measured at some optical frequency A. 

Having the singular values of the Jones matrix [J] or eigenvalues of [H] = [J]t[J] 
matrix the transmission loss can be computed as follows [80]. 

lout = (E 0 ut\E 0ut ) = (^j n |[J] t [J]|^ irl ) = 

= <# in |A|# in ) = ^^I in (4.8) 

Which can be re-fornrulated in the more frequently used dB scale as: 

hut mi P1 + P2 

— = 101og 10 —-—. (4.9) 

J-in ^ 

Similarly, the polarization-dependent loss, which is the magnitude of the difference 
between the maximum and minimum device transmission over all possible input po¬ 
larization actually corresponds to the difference in the loss measured along the system 
principal axes U\ and u 2 . It is defined as follows: 

PEL = 10| log 10 — | (4.10) 

P2 
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Alternatively we can introduce a parameter called degree of polarization [3]: 

P = Pl — P2 (4.11) 

Pi + P 2 

The eigenvalues pk = can be interpreted as an observable transmission loss for 
the case when incident to an optical system electric vector E is align with the system 
principal axes, i.e. E || Uk- The PDL parameter accounts for the difference in the 
loss measured along the system principal axes Ui and u 2 . 


eigenvalues pi 'll Pi‘Pi eigenvalues Pi ■ Pi 





Figure 4.7: Eigenvalues pi and p 2 and the angle 0 between the geometrical axes of 
the system x, y and the coordinate system defined by the principal axes ui, u 2 , before 
a) and after b) the eigenvalues were reordered. (The data were obtained by means of 
the OVA 5000 Luna Technologies.) 


In addition to the eigenvalues of the system transmission matrix (plotted as a 
function of wavelength in Figure 4.7), we can also find the angle of rotation (0) of the 
principal axes about the device optical axis (z) and plot it as a function of wavelength 
(bottom frame of Figure 4.7). From the mathematical point of view, both eigenvalues 
are roots of a quadric polynomial and can be ordered arbitrary, usually in descending 
order: this is why p\ is always larger than p 2 in Figure 4.7a, hence the eigenstates are 
interchanged when they cross. Because of this effect, the corresponding eigenvectors 
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are interchanged as well, therefore the angle 0 experiences 7t/2 shifts at the crossing 
points, as shown on the bottom panel of Figure 4.7a. Therefore, a strategy for 
restoring the individual transmission spectra along the principal axes is to interchange 
the eigenvalues and principal axes every time the 0 = 7t/2 jumps are detected. 

The result of such reordering is shown Figure 4.7b. Now the transmission spectrum 
corresponds to linearly polarized light with its electric held vector aligned with each 
principal axis of the system and the rotation angle jumps are eliminated. 

In general, the principal axes of an optical system are not necessarily fixed with re¬ 
spect to a reference frame but may depend on the wavelength, as shown in Figure 4.8 
for a TFBG sensor. Indeed, by changing the optical frequency it is expected that an 
optical system would operate differently. The global behavior of the principal axes, 
along the whole operational range of the TFBG sensor, is shown in Figure 4.8b, and 
the corresponding insertion loss in Figure 4.8a. To remove the noise only the points at 
which the difference between eigenvalues is noticeable (|p 2 — pi\ > 0.05) are plotted. 
The noise arises from the fact that at the points of eigenvalues crossing the trans¬ 
mission matrix becomes degenerate, has the eigenvectors are not well defined, and 
oscillate rapidly with respect to the geometrical axis, as can be seen from Figure 4.8b 
at the points of crossing. 

Figure 4.8b shows that the system possesses significant polarization asymmetry, 
or birefringence, in the wavelength range A G [1545 — 1575]. Although the principal 
axes are globally stable, near 75° degrees relative to the reference frame of the LUNA 
interrogation system, locally they experience small oscillations, of about 8° degrees 
about the optical axis. These oscillations are related to the resonances observed in 
the spectrum (Figure 4.8a), and reflect a wavelength dependent birefringence. 

We also note the observed alternation between the peaks (Figure 4.8a), which 
can be explained by the fact that the alternating peaks have different azimuthal 
symmetries and polarizations [9]. 
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Figure 4.8: Rotation of the principal axes of TFBG device as a function of optical 
wavelength. (The data was obtained by means of the OVA 5000 Luna Technologies.) 


4.2.2 Measurements along geometrical axes of an optical 
system 

Instead of choosing principal axes as a reference coordinate system, geometrical 
axes can be chosen alternatively. In the case of TFBG it is convenient to choose 
orthogonal axes such that one of the axis is normal to the fibre axis and lay on the 
plane parallel to the grating blades, as shown in Figure 3.36. Usually the spectra 
along geometrical axes are obtained by introducing an external linear polarizer, such 
that state E in can be align along the required axis. In Figure 4.1 such orthogonal 
states along x and y axes are shown on the density plot and denoted as I x and I y , 
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respectively. 


4> = 0° (aligned axes) 



Wavelength, (n m) 



Wavelength, (n m) 


Figure 4.9: Transmission loss along the TFBG system principal axes and its geomet¬ 
rical axes for perfectly align coordinate systems 0 = 0° and rotated by 45° degrees 
(j) = 45°. The intersession loss I x>y and eigenvalues pi .2 are given in linear scale. 


The difference in transmission loss along the principal axis of the TFBG sensor and 
its geometrical axes is shown in Figure 4.9. ft can be clearly seen that in the case when 
geometrical axes are align with the principal axes of an optical system both approaches 
yield almost identical result. The small difference comes from the aforementioned 
wavelength dependent oscillation of the principal axes (which cannot be compensated 
for in the direct measurement), and to a possible slight change in polarization state 
between polarizer and the grating (since a non-polarization maintaining fibre is used). 
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The following section will demonstrate that in spite of this small inaccuracy, the 
parameters extracted from the Jones matrix data provide excellent spectral sensitivity 
results for refractometric sensing. 

Nevertheless, although the singular value decomposition, or egen decomposition, 
provide an elegant way to introduce coordinate axes along which the system can be 
studied, it might be desirable to fix the reference axis permanently to an optical 
system geometry, for example in the cases of TFBG sensor, when it is known that 
the system behaves physically different along particular axes. 

4.2.3 Extracting transmission spectra measured along 

geometrical axes from the Jones matrix or the Stokes 
vector data 

In this section we describe an approach allowing to extract transmission loss spec¬ 
tra along the given geometrical axes from the Jones matrix or the Stokes vector 
without introduction of external polarizer. Thus a single measurement by OVA can 
provide a complete information about transmission loss along all possible geometrical 
axes. 

The Jones matrix analysis 

Transmission spectrum of linear polarized light, with the electric field vector E 
align along a given geometrical axis can be extracted directly from the Jones matrix, 
measured for example by means of OVA 5000 from Luna Technologies. 

The instrument was designed to provide a full characterization of the polarization- 
dependent transmission properties of optical fibre devices as a function of wavelength 
in the form of the Jones matrix elements. The spectral accuracy of the OVA is 
+/ — 1.5pm and its insertion loss accuracy is +/ — 0.1 dB, over a wavelength range 
from 1525 to 1610 nm. 

The OVA 5000 Luna is capable of capturing four complex transfer functions (am¬ 
plitude and phase) of a fibre as a function of wavelength and polarization. The 
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incident and transmitted light can be represented in terms of the electric field vector, 

, known as a Jones vector [79], where the field 

components E x , E y are projections of the electric held vector E on the x, y axes. The 
optical TFBG sensor then can be represented as a four port device with two input 
and output states of polarization. These four transfer functions can be assembled 
into a Jones matrix [79], giving the complete characterization of a device under test. 

Eout — [J]Ein (4-12) 


written as a column vector: E = 


E x 
. E v 


[■/(A)] 


a(A) b{ A)\ 
c(A) d(X)J 


(4.13) 


Here, J(A) is the Jones matrix consisting from the four complex transfer functions 
a(A), 6(A), c(A), d(X) of the optical frequency A. 

The Jones matrix of an optical device under test connected to a linear polarizer 
is give by: 

[d dev+pol] ['J pol ] [ '7 dev ] • (4-14) 


Here [Jdev] is the measured by an optical vector analyzer Jones matrix of a bare 
device, and [J po i] is the Jones matrix of linear polarizer known from theory. 

The Jones matrix of a linear horizontal polarizer [Jdev(9)\, rotated by 9 angle with 
respect to the optical device axes, is given by: 


[Jvoim 


im) 


1 0 
0 0 


[R(-9)] = 


cos 2 (6>) cos(0) sin(0) 
sin( 6 l ) 008 ( 6 *) sin 2 (6*) 


where [77(0)] is the rotation matrix: 


im) 


cos (9) —sin (9) 
sin(6*) cos(6*) 


(4.15) 


(4.16) 
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rotating the coordinates system so that the horizontal polarizing element has the 
simplest representation and then rotating the system back down into the original 
system of coordinates. 

Tims, on optical device with a known Jones matrix [Jdev\, connected to a linear 
polarizer rotated by 9 angle, is described by the Jones matrix: 


[ Jdev+pol ( @ ) ] 


cos 2 (9) 
sin( 6 *) cos( 0 ) 


008 ( 6 *) sin( 6 *) 
sin 2 (9) 


[J< 


dev I • 


(4.17) 


Finally the transmission spectrum 1(9, 9), for a given angle 9 of the linear polarizer, 
can be computed as 


/(M) 


10 log 10 


Pi(9, A) + p 2 (9, A) 


(4.18) 


where p x (9, A) and p 2 (9, A) are the eigenvalues of 77(0, A) = [J de v+ P oi(Q, A )] j [J d ev+ P oi(0, A)] 
matrix. 


The Stokes vector analysis 


In a similar way we can use the Stokes vector. The Stokes vector can be measured 
with a polarizer controller or a special instrument, such as JDS Uniphase SWS-OMNI- 
2 system. 

A beam of light can be completely described by the four parameters [81, 82], 
represented in the form of the Stokes vector: 




( 7(0°)+ 7(90°) ^ 

Si 


7(0°) -7(90°) 

s 2 


7(45°) -7(135°) 

W 


\ Irhs ~ Ilhs ) 


(4.19) 


here I(9) is the intensity of light polarized in the direction defined by the angle 9 in 
the plane perpendicular to the direction of light propagation, and Irhs, Ilhc are the 
intensities of right- and left-handed polarized light, respectively. 

Since the Stokes parameters are dependent upon the choice of axes, they can be 
transformed into a different coordinate system with a rotation matrix [ R]. Considering 
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that a second coordinate system is obtained by rotating the original coordinate system 
about the direction of light propagation on the angle 9 , we can write [83] 


S'(e) = [fl(9)]s, 


(s' \ 

‘-’O 


10 00 





0 cos(2 9) sin(20) 0 



S' 2 


0 — sin(2 9) 008 ( 20 ) 0 



UJ 


0 0 0 1 


w 


S„ \ 


Si cos(20) + S 2 sin(20) 

— S\ sin(2 9) + S 2 cos(2 9) 

V s » J 

here So, Si, S 2 and S 3 are the Stokes parameters. 

Now considering that 


(4.20) 


(4.21) 


7(0°) +7(90°) = S 0 , 

7(0°)-7(90°) = Si cos(20) + S 2 sin(20), (4.22) 

we conclude that the transmission loss spectra along the two orthogonal axis, rotated 
by the angle 9 with respect to the original system of measurements, are given by the 
following expressions: 

7 x (9, A) = i(5 0 (A) + -S'i.(A) cos(20) + S 2 (A) sin(20)), 

I y (9 , A) = ^(S 0 (A) - S'i(A) cos(20) - S 2 (A) sin(20)). 

(4.23) 


Hence, if the Stokes vector is known in one coordinate system, the transmission spec¬ 
trum of linearly polarized light with the electric field E align along a given geometrical 
axis, defined by the angle 9 , can be recovered. 
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We should also note that the Jones matrix and Mueller matrix representations are 
connected, and can be expressed as [84]: 


M = U{J®J)U~\ 


(4.24) 


where J ® J is the direct product of Jones matrices, and matrix U is given by 


U 


10 0 1 

10 0-1 
0 110 
0 i —i 0 


(4.25) 


Although the two approaches are related the phase information is available only in 
the case when the Jones matrix is known. The phase information allows to determine 
additional phase related parameters, such as Group Delay (GD). 


4.3 Polarization-based detection of small 

refractive index changes with TFBG sensors 

In this section we investigate which polarization-based measurement techniques 
provide the best signal-to-noise ratio when TFBG sensors are used to detect small re¬ 
fractive index changes, and use the special case of a TFBG coated with gold nanorods 
(as described in [3]). This choice is made because the polarization dependence of 
waveguide-type sensors is much enhanced when metal coatings are used. As indi¬ 
cated earlier, the interaction of guided waves with metal interfaces depends strongly 
on whether the electric fields of the waves are tangential or normal to the metal bound¬ 
ary. It was further mentioned that when the core-guided input light of a TFBG is 
linearly polarized along the principal axes, the electric fields of high order cladding 
modes are either tangential or radial (hence normal) to the cladding boundary. To be 
precise, y-polarized input light (corresponding to light polarized parallel to the plane 
of incidence on the tilted grating fringes, as shown in Figure 3.36, i.e. P-polarized) 
couples to radially polarized cladding modes while x-polarized light (perpendicular 
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to the tilt plane, or S-polarized) couples to azimuthally polarized cladding modes 
(tangential to the boundary). Further polarization effects arise with non-uniform 
metal coatings, since they have boundaries that are both tangential and radially ori¬ 
ented relative to the cylindrical geometry of the fibre [85]. It is therefore desirable 
to carry out two transmission measurements along the principal axis to observe di¬ 
rectly such polarization effects on the strengths and positions of the cladding mode 
resonances [ 86 ]. On the other hand, it has been shown here that a Jones matrix 
measurement can provide this information as well, in addition to other parameters 
of interest, such as PDL. While the PDL spectrum “hides” the physical effects re¬ 
sponsible for difference in transmission due to the different polarization states, it has 
been shown in the past to yield excellent limits of detection for surface plasmon res¬ 
onance based TFBG sensors [15]. We now proceed to compare the signal noise for 
refractive index measurements by the various polarization-dependent data extraction 
techniques. 

As shown in Figure 4.10 for a TFBG coated with a sparse layer of gold nanorods 
and immersed in water, the PDL parameter provides the absolute value of the dif¬ 
ference between resonances observed in the transmission spectra measured along the 
principal axes. The relative position of the peaks, their amplitude, and their width 
become convoluted in the PDL parameter, which provides instead a maximum lo¬ 
cated somewhat in between the individual resonance maxima and a zero on either 
side corresponding to the wavelengths where the spectra cross each other. 

When the refractive index of the medium surrounding such TFBG changes the 
waveguiding characteristics of the cladding are modified and the resonances observed 
in the transmission spectrum change accordingly. Therefore, to detect changes in 
refractive index we can either follow the amplitudes and positions of individual reso¬ 
nances [ 86 ] or of the PDL features. Here, the sensor was immersed in water and the 
refractive index was incrementally increased in steps of An = 1.517 x 10~ 4 by adding 
10 /il of ethylene glycol ( C 2//4 ( 0 / 7)2 ) to 5 ml to the water. The impact of each 
increase in refractive index on all parameters of interest is shown in Figure 4.11 for 
a typical slice of the spectrum (it was noted in [3] that there was little difference in 
wavelength shifts across the TFBG spectrum for this device). 
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Pi <P2 , (dB) 




Wavelength, X (nm) 

Figure 4.10: Two eigenvalue spectra (individual transmission along the principal axes) 
and corresponding polarization-dependent loss (PDL) parameter. The peak in PDL 
spectrum is denoted by “1” and zeros by “2”. Here the more common dB scale is 
used. 


The refractive index change was chosen to have a relatively small value of ( An = 
1.517 x 1(T 4 ) to test the sensor detection limits and a linear fitting was used because 
the sensor response for such small changes is expected to be linear. The standard 
deviation of the errors from the linear fit was calculated in the usual manner by: 


SD = 


\ 


N 




(4.26) 


i —1 


here // is the expected value of x, and x = y appr — Vdata is the difference between the 
measured data ijdata and its linear approximation y appr ■ 




















Chapter 4. Experimental polarization-based sensing with TFBG 


129 


Mnm) 'x peak ' SD= 2.83 pm (nm) l y peak, SD = 2.92 pm 




Figure 4.11: Position of the resonance a) in the transmission spectra I x and I y of 
light polarized along x and y geometrical axes here aligned with the principal axes, 
and b) in PDL spectra (the maximum and zero values of PDL are detected), as a 
function of refractive index change. The continuous lines represent the least square 
approximation to the measured data, and SD is the standard deviation from the 
linear approximation 

The results of the fits show that the most accurate detection of the refractive index 
change is achieved with the zeros in the PDL spectrum, as this measurement provides 
the smallest standard deviation of SD = 1.40 pm. The worst result were obtained 
with the PDL peak (SD = 6.44 pm), while the detection of individual polarized res¬ 
onances provides an intermediate value of the standard deviation (essentially equal 
to 3 pm), ft is not surprising that zeros of PDL should provide the most accurate 
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results as they consist essentially of a differential measurement between two spectra 
with well-defined crossing points that occur on the sides of the individual resonances, 
where the spectral slope is highest. On the other hand, at the expense of an in¬ 
crease in noise by a factor of approximately 2, other effects such as differences in the 
change in the resonance amplitude for the two polarizations (which can be linked to 
differential loss or scattering) can be studied when the principal axis spectra are used. 

4.4 Conclusion 

In this Chapter we investigated the use of Jones matrix and Stokes vector based 
techniques and polarization analysis to extract information from optical fibre sensors 
in non-polarization maintaining fibres. In particular we showed how to calculate 
transmission spectra with electric holds E aligned with the system principal axes or 
along any other system axes. We also proposed the method of extraction spectra of 
light polarization along a given geometrical axis either from the Jones Matrix or the 
Stokes vector data. In the case of a TFBG inscribed in a non-polarization maintaining 
fibre for instance, the principal axes of the system are determined by the direction of 
the tilt of the grating planes (and its perpendicular). The transmission spectra for 
light polarized in the tilt plane (P-polarized) or out of the tilt plane (S-polarized) can 
thus be extracted without having to separately align a linear polarizer upstream from 
the grating and hoping that the polarization remains linear between the polarizer 
and the TFBG. Polarization-resolved spectra can also be obtained at a faster rate, 
for applications in chemical deposition process monitoring for instance, without the 
need to line up or rotate a polarizer between each measurement. 

The transmission spectra of P- and S-polarized light were compared with the 
TFBG spectral response along the principal axes. A small oscillation of about 8 de¬ 
grees in the orientation of the principal axes as a function of wavelength was observed. 

Furthermore, it was determined that the best TFBG sensor detection limits are 
achieved when zeros in the PDL spectrum are followed, mainly because of the sharp¬ 
ness of crossing points occurring on the steep sides of the individual resonances, but at 
the expense of the direct observation of the sensor response to polarized light (which 
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are also available from the Jones matrix data). In the latter case, while the standard 
deviation of the spectral sensitivity is doubled relative to PDL measurements, addi¬ 
tional information becomes available regarding the influence of the measured medium 
on cladding mode loss, allowing further uses from the TFBG data sets [87]. 



Chapter 5 


Optical properties of materials 

In the presented work we investigate methods of increasing TFBG sensor sensi¬ 
tivity by coating its surface with various types of coatings, such as metal films and 
films consisting of nanoparticles. 

Before running numerical simulations, the optical properties of materials should 
be defined. In this chapter we will focus our attention on optical properties of metals 
since it is known the that sensitivity of an optical sensor can be increased by excitation 
of surface plasmon resonances (SPR) observed in metal films and nanoparticles of a 
proper geometrical dimension [88, 89, 90, 91]. We also consider various surface mor¬ 
phologies and review methods that allow us to account for the dipole-dipole radiative 
interaction between elements of a him or between closely deposited nanoparticles. 

5.1 The methods of measurement of optical 
constants 

It is not a trivial task to choose the correct optical parameters for thin metal films 
or nanoparticles. Various sources of information do not always provide consistent 
results. The theoretically computed permittivity of metals is not always consistent 
with the experimental measurements, and experimental results sometimes are limited 
by a particular method. 
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In this section we provide a brief review of various methods available for measuring 
of optical constants and discuss limitations imposed by a particular method, so that 
we can choose a reliable data source for further computation. Next we will consider 
models for mixed media needed to describe rough films resulting from nanoparticle- 
based coatings. 

The optical constants n and k of a medium at frequency uj can be determined either 
by (1) measuring n and k independently at a given oj or (2) one of the parameters 
(n or k) should be measured over the entire frequency range, then the remaining 
parameter can be deduced from the measured data. 

Several alternative methods were proposed as well, including the Drude ellipsom- 
etry method introduced in 1889 [92] in which the relative amplitude and phase shift 
are measured simultaneously, however, the method is significantly affected by the 
sample surface quality [93]. 

5.1.1 The methods based on single parameter 
measurements 

Historically, this approach was based on Bode’s electric network theory introduced 
in 1945, where it had been shown that the attenuation and phase at the output of 
an electric network are not independent of one another. In Robinson’s paper from 
1952 [94] the absorption spectrum was derived from a normal incidence reflection 
spectrum. The phase change of a reflected wave was deduced from the curve of the 
reflection attenuation as a function of frequency, and the optical constants n and k 
were then determined with the aid of a Smith chart. 

A modern treatment of the problem can be found in [95, 96, 97]. The usual optical 
procedure is to measure the reflectance at normal incidence over as wide range as 
possible, and then use, for example, the free-electron extrapolation for the reflectance 
outside the range. Once the reflectance is approximated over the range from zero to 
infinity the Kramers-Kronig (KK) integration might be implemented: 
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X2 M 


XiM 



(5.1) 


Here the integration is done in the sense of the Cauchy principal value, defined as the 


P f and x(cu) is a complex analytical function, vanishing faster than j^- as |co| —* oo, 


XiM := Re[x(u)\ and y 2 (w) := Im\x(u)\. 

Thus the optical constants n and k can be connected through the Kramers-Kronig 
relations: 


x(cj) = n{uj) = n{uj) + ik(u) = 
- xiM+*»(«) 


(5.2) 


However, the results of such calculations are often quite sensitive to the shape of 
the extrapolated spectrum, therefore the extrapolation of the reflectance at both ends 
of the spectrum is crucial. While a reasonable low-energy extrapolation is possible by 
simply assuming that R = 100%, additional experimental information is needed to 
make a reasonable high-energy extrapolation. The measurements have to be extended 
further into the vacuum ultraviolet region in order to determine the n and k in the 
visible band. The necessity of such measurements is the major disadvantage of this 
method [97, 98]. 

5.1.2 The methods based on simultaneous measurement of 
both parameters 

Simultaneous measurement of both parameters is usually difficult in practice, es¬ 
pecially in the region of strong absorption. To overcome this problem several new 
techniques were proposed and successfully implemented in [98, 4], Instead of a single 
normal incidence reflectivity measurement, the parameters for oblique incidence at 
several angles and various polarizations were measured. The right choice of the inci¬ 
dent angle and polarization state can significantly improve the results in the region 
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where normal incidence measurements are inaccurate [98]. The proposed method is 
also less sensitive to the surface roughness, unlike the ellipsometry method. 

First the reflectance R and transmittance T of a thin film deposited on a transpar¬ 
ent substrate were calculated for the normal and oblique incidence, for s and p polar¬ 
ization, and different him thicknesses. The calculation involves solution of Maxwell’s 
equation boundary-value problem from which the reflectance R and transmittance T 
can be expressed in terms of the material constants n and k. 

Next the obtained function R(k,n ) and T(k,n ) were inverted to obtain n and 
k expressed in terms of the measured R and T. Since the optical constants n and 
k are independent of the boundary conditions, the n and k values were expected to 
depend only on the incident light wavelength (which excites the electronic transition), 
but not on polarization and angle of incidence (if material is isotropic) or on the 
him thickness. Therefore at each wavelength, in the range of interest, reflectance R 
and transmittance T at specihed incident angles and polarizations were measured. 
Then the him thickness was measured by means that were independent of the R 
and T measurements and compared with the result deduced directly from R and T 
measurements [98]. The process of inversion requires solution of several nonlinear 
equations, and a number of iterations till the convergence is achieved [98]. 

5.2 Optical constants of metals 

In our calculation we are going to use optical constants obtained in [4] by the 
method described in the previous section [98], where the optical constants n and 
k were obtained for copper, silver, and gold from reflection and transmission mea¬ 
surements. The films were created by vacuum-evaporation with hlm-thickness in the 
range of 185 — 500 A, and the measurements were conducted in the spectral range of 
0.5 - 6.5 eV. 

It was observed that the results were independent of him thickness only above 
a certain critical thickness of about 250 A, and were unchanged after vacuum an¬ 
nealing or aging in air. The optical properties of evaporated thin films have been 
found to be the same as for bulk materials for the thickness of the films greater than 
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about 300 A [4], 


5.2.1 Free electron approximation, the Drude-Sommerfeld 
model 

The obtained experimental data can be approximated with the free-electron model. 
Assuming that the motion of electrons is confined to a region much smaller than the 
wavelength we can implement the Drude-Sommerfeld model [99, 92], The model does 
not include a restoring force, assuming free electrons, thus an equation of motion of 
a single electron can be written in the following form: 

m 0 ■ d 2 t x(t) + 7 ?vy • d t x(t) = F ext {t ), (5.3) 

where m 0 is the effective optical mass of electron, 7 = - is the macroscopic damping 
constant due to the dispersion of the electrons caused by the crystalline structure, 
r is the relaxation time and F ext (t ) is the external force applied to the electron. 

Assuming harmonic excitation: F ext (u) = — eE(uj ) and taking the Laplace trans¬ 
form of equation (5.3) we have: 


x{oj) = 


-E{u) 


m 0 cu 2 + iuj'y 

The macroscopic polarization can be written in the following form: 


(5.4) 


Ne 2 1 

P(u) = N ■ p(u) = —Ne ■ x{u) =--- E(u) = 

m 0 ur + i^/oj 

= e a x(uj)E(u) = 

= e a (e(ui) — 1) E(u>), (5.5) 

where: x is die electric susceptibility, N is the number of conducting electrons per unit 
volume (density), eo is the electric permittivity of free space and r is the relaxation 
time. 

Excluding E(uj) from equation (5.5), the expression of e(cu) can be obtained in 
terms of effective optical mass and macroscopic damping of a free electron. We denote 
it as 
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sf(u>) = 1 




u 2 + i'yu 


here c u p = 


yj is the plasma frequency. 


The complex dielectric constant 


(5.6) 


e — ei + u 2 


(5.7) 


and the complex index of refraction, (measured in experiments) 


n — n + ik 


are connected by the relation e = n 2 , so that 


(5.8) 


2 / 2 
ei = n — k , 

62 = 2 nk. 


(5.9) 


Considering (5.6) it is possible to separate e-f into its real and imaginary parts: 


e{ = 1 


J - 

to - 


uj 2 t 2 


1 + urr 


2~2 ’ 




(5.10) 


cn(l + uj 2 t 2 ) 

We note that free-electron approximation to the dielectric function is determined 
by the relaxation time r and the the plasma frequency u p defined by the electron 
optical mass m 0 . 


5.2.2 The near infrared band 

For metals at near-infrared frequencies we can assume uj y. Considering equa¬ 
tions (5.6) and (5.10) we can write: 

2 = ,-_> 2 
,2 


£ /. i_s = 

1 -2 A 2 ’ 


, uj 2 
f p 

e 2-5“ 

l jJ a T 


W 

A 2 ’ 

p 


(5.11) 
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where 


1 _ Ne 2 

A 2 7T77l 0 C 2 ’ 

T = 27 TCT. 


(5.12) 


The values of n(oj) and k(oj) can be measured experimentally, hence the values of 
and e 2 can be calculated using (5.9). Assuming that the experimental results at 
near-infrared bands can be satisfactorily approximated with the free-electron model 
we set 

e = e / . (5.13) 

Then the optical mass m 0 can be determined from the experimental results for e\ 
from the slope of a plot of —e vs A 2 . Next, using the expression for e 2 and plotting 
e 2 /A vs A 2 we can determine r from the slope of the graph in the infrared band. Such 
derivation was conducted in [4] where the effective optical mass and relaxation time 
were obtained, ft was shown that the effective optical mass and relaxation time can 
be considered to be constant for the whole near-infrared band. 

The results for copper, silver, and gold are shown in Table 5.1 [4]. The results for 
m 0 are relatively accurate since in the infrared k^> n and can be measured precisely. 
The error in the r measurements is larger due to the error in the n measurements in 
the infrared band [4]. 

Table 5.1: Optical masses and the relaxation times for copper, silver and gold. [4] 


Metal 

mo (electron masses) 

r x 10 15 (sec) 

Copper 

1.49 ±0.06 

6.9 ±0.7 

Silver 

0.96 ±0.04 

31 ± 12 

Gold 

0.99 ±0.04 

9.3 ±0.9 


We can conclude that the properties of metals at the near-infrared band can be 
approximated with the free-electron model, defined by the parameters m 0 and r, 
which can be obtained from the experimental measurements of the optical constants 
n and k. 







Chapter 5. Optical properties of materials 


139 


5.2.3 The visible band. The interband absorption 

In the visible and near-ultraviolet regions Drude’s free-electron theory fails and 
the interband absorption should be taken into account. The absorption in the visi¬ 
ble and ultraviolet band is significantly influenced by the transitions from the com¬ 
pletely occupied d bands to an empty state above the Fermi level in the conduction 
band. Moreover, these transitions depend on the nanoparticle size. The interband 
absorption peak becomes sharp and its peak position shifts towards the lower energy 
side [ 100 ]. 

We are going to review the size-induces changes in the d band in the following 
chapter, here we simply state the significance of the interband absorption. The in¬ 
terband contribution to the imaginary part of the dielectric constant can be obtained 
by subtracting the free-electron contribution value el from the experimentally deter¬ 
mined value of 62 - The comparison between the experimental data and theoretical 
prediction, based on the classical free-electron Drude theory, was conducted in [4], 

The interband contribution was also compared with the theoretical prediction 
based on the band-structure model, where the joint density of states and the transition- 
probability matrix elements throughout the Brillouin zone were taken into account [101, 
102, 103], 

The discussion of optical constants of Ag and Cu in terms of free-electron ef¬ 
fects, interband transitions and collective oscillations can be found in Ehrenreich and 
Philipp [104], 

We conclude that the free-electron expression of e is useful only for photon en¬ 
ergies below a threshold energy and can be applied for the near-infrared band only. 
Above this threshold energy the form of the 62 curve depends on the specific material 
band structure, and also exhibits size-dependent properties [ 100 ]. 

Assuming that interband contribution in known either from experiments or theory 
we can write 

e(cj) = e in tra(u) + e f (u) (5.14) 

where e(oj) is the dielectric function of a metal, e^(a;) is the free-electron contribution 
and €intra(w) is the interband contribution. 
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Although the interband contribution can be satisfactory calculated only in the 
framework of quantum mechanics, or obtained from the experimental measurements, 
the classical approximation is still possible. The electrons in the d band are not free 
as valence band electrons, but rather bounded by the lattice ions. Correcting the 
classical free electron model (5.3), by introducing a linear restoring force we come 
to the Lorentz-Drude oscillator model, based on the damped harmonic oscillator 
approximation: 


m d ■ d 2 t x(t) + m d 7 • d t x{t) + m d u>}x(t) = F ext (t), (5.15) 

thus the single induced dipole moment is: 

e 2 1 

p(u) = exito) =-^---:— E(oj). (5.16) 

v ’ v ’ m -u 2 d + u 2 + zcuy v K ’ 

The macroscopic polarization P(u>) depends on the model describing dipole-dipole 

interaction. The discussion on this subject can be found in [105]. 

In practice it is not always possible to approximate interband absorption with 
one specific type of electron (which has an effective mass m d and resonates at some 
eigenfrequency c o d ). Usually several different types of electrons with different effective 
masses bounded by different restoring forces should be taken into consideration. 

Thus for each resonance in a given absorption spectrum the restoring forces and 
the fraction of each type of electrons can be chosen for accurate approximation. Such 
investigation was conducted in [5] were the six-oscillator Lorentz-Drude model was 
used to fit optical functions of eleven widely-used metals in optoelectronic. The 
parameters were chosen for the best fit of experimental data, obtained by the method 
described in the previous section [4], 

In accordance with Lorentz-Drude oscillator model, a complex dielectric function 
e(cu) can be expressed [106, 5]: 

e(u) = e f (u) + e b (u), (5.17) 


where the free-electron or Drude model is represented by the first term: 


e f (u) 


= 1 - 


fo^l 


Ld 2 — icuTn 


(5.18) 
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and the interband part of the dielectric function is represented by the bounded Lorentz 
electron model, similar to the model used for insulators: 


N 


M = -£ 


f 


3 =1 






iouTj 


(5.19) 


Here u p is the plasma frequency, N is the number of oscillators with resonant fre¬ 
quency u)j, and oscillator strength fj, and Tj = d- is the macroscopic damping con¬ 
stant of the jth oscillator (To corresponds to the plasma damping constant ). 

The obtained results can be summarized in Table 5.2 [5] (where c Oj is given in 
electron volts, and Tj in sec" 1 ). 
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Table 5.2: Values of the Lorentz-Drude Model Parameters [5]. 



Ag 

Au 

Cu 

A1 

Be 

Cr 

Ni 

Pd 

Pt 

Ti 

w 

(dp 

9.01 

9.03 

10.83 

14.98 

18.51 

10.75 

15.92 

9.72 

9.59 

7.29 

13.22 

Ui 

0.816 

0.415 

0.291 

0.162 

0.100 

0.121 

0.174 

0.336 

0.780 

0.777 

1.004 

i02 

4.481 

0.830 

2.957 

1.544 

1.032 

0.543 

0.582 

0.501 

1.314 

1.545 

1.917 

U) 3 

8.185 

2.969 

5.300 

1.808 

3.183 

1.970 

1.597 

1.659 

3.141 

2.509 

3.580 


9.083 

4.304 

11.18 

3.473 

4.604 

8.775 

6.089 

5.715 

9.249 

19.43 

7.498 

^5 

20.29 

13.32 

0 

0 

0 

0 

0 

0 

0 

0 

0 

r G 

0.048 

0.053 

0.030 

0.047 

0.035 

0.047 

0.048 

0.008 

0.080 

0.082 

0.064 

r, 

3.886 

0.241 

0.378 

0.333 

1.664 

3.175 

4.511 

2.950 

0.517 

2.276 

0.530 

r 2 

0.452 

0.345 

1.056 

0.312 

3.395 

1.305 

1.334 

0.555 

1.838 

2.518 

1.281 

r 3 

0.065 

0.870 

3.213 

1.351 

4.454 

2.676 

2.178 

4.621 

3.668 

1.663 

3.332 

r 4 

0.916 

2.494 

4.305 

3.382 

1.802 

1.335 

6.292 

3.236 

8.517 

1.762 

5.836 

r 5 

2.419 

2.214 

0 

0 

0 

0 

0 

0 

0 

0 

0 

fo 

0.845 

0.760 

0.575 

0.523 

0.084 

0.168 

0.096 

0.330 

0.333 

0.148 

0.206 

fl 

0.065 

0.024 

0.061 

0.227 

0.031 

0.151 

0.100 

0.649 

0.191 

0.899 

0.054 

h 

0.124 

0.010 

0.104 

0.050 

0.140 

0.150 

0.135 

0.121 

0.659 

0.393 

0.166 

h 

0.011 

0.071 

0.723 

0.166 

0.530 

1.149 

0.106 

0.638 

0.547 

0.187 

0.706 

U 

0.840 

0.601 

0.638 

0.030 

0.130 

0.825 

0.729 

0.453 

3.576 

0.001 

2.590 


5.646 

4.384 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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5.2.4 Dispersion curves 

In this section, as an example, we plot the real and imaginary parts of the refractive 
index n(cu) = n(co) + ik(co) of Au, Ag, A1 and Cu metals. The optical properties of 
metals were described phenomenologically using the Lorentz-Drude model, in which 


parameters of six oscillators were fitted for 
The results are shown in Figure 5.1. 



Figure 5.1: Optical constants n and k of 
energy. 


the best consistency with experiments 



Frequency © (eV) 

, Ag, Cu and A1 as a function of photon 


It can be seen from Figure 5.1 that the simple free electron Drude model can be 
applied only below the visible band for Ag and Au, where is in the case of A1 and Cu 
the interband absorption resonances are present even in the IR band. 

Once the dispersion curves of various metals are known, we can proceed with 
simulation of various metal-based coatings. 
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5.3 Optical properties of mixtures and rough 
surfaces 

In this section we review a simple yet efficient approach allowing us to describe a 
nonuniform metal coating, rough surfaces and surfaces coated with nanoparticles. 

5.3.1 Local field effects and effective medium theory 

Probably the simplest way to calculate dielectric permittivity of a given material 
with inclusions of another material is to apply the effective medium theory. The 
effective medium theory describes macroscopic properties of a medium based on the 
properties and the relative fractions of its components. 

In 1909 Lorentz [107] had shown that the dielectric properties of a substance can 
be related to the polarizabilities a with the Clausius-Mossotti (or Lorentz-Lorenz) 
relation [108, 109, 40]: 


3 e — 1 . 

a = — -, (5.20) 

lVe + 2 v ’ 

where N is the number of dipole particles per unit volume and e is the dielectric 

permittivity. The derivation goes as follows [110, 40]: an external field E ext induces 

local molecular dipole moments pj at each jth site, proportional to the local field 

Ei oca i , with the proportionality constant a - called the atomic polarizability: 


P OiEiocai. 


(5.21) 


In the general case cc is a complex number, meaning that the polarization may be 
shifted in phase with the external electric field, dependent on the external field fre¬ 
quency. The polarization can be a tensor for a non isotropic material. The local held 
Ei OC ai is not equivalent to the external held E ext and should be corrected 1) by taking 
into account the held from all other dipoles (excluding the one under consideration) 
and 2) in some cases by including the dipole radiation, if the external held oscillates 
at a sufficiently high frequency. 
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The first effect is considered by subtracting a microscopic sphere from the con¬ 
tinuous medium and finding the electric field inside the formed imaginary cavity. 
The sphere diameter has to be chosen much smaller than the wavelength, hence the 


average field may be assumed to be spatially homogeneous. The resulting field is 


the local field Ei oca i (which corrects the E ext field) and is directly connected to the 
induced local dipole moment p. The macroscopic polarization P is simply a sum of 
all the local microscopic dipole moments. Then macroscopic observables, such as the 
dielectric function of a medium, can be easily deduced from the known connection 
between an external field E ext and the correspondent macroscopic polarization P. 

The second effect (the dipole radiated field) can be neglected in many cases, 
and the result can be viewed as a zero-frequency limit [40]. If the effect of the 
dipole radiative interaction is significant it can be accounted with the discrete dipole 
approximation method. 

In the presented work we are interested in the calculation of the dielectric per¬ 
mittivity of particles made of the same material (the inclusions) which are embedded 
in another material (the host). The derivation of the general optical mixing formula 
can be found in [105, 40, 111] and written in the following form: 


e- e h 



(5.22) 


e/i + (e — th)L 


where 

e - is the effective dielectric function. 

6j - the dielectric functions of the constituents. 

6h - is the dielectric function of the host material. 

L - is the depolarisation factor, defined by the morphology. For spherical inclusions 
L = 1/3. 

Pj = Vj/V - is the filling factor of the constituents, also called the volume fraction. 

V - is the full volume occupied by mixtures, 

Vj - is the volume fraction occupied by jth constitute. 

By making several assumptions the general form of equation (5.22) can be reduced 
to a set of particular cases, best suited for certain types of media, such as aerosols, 
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porous media and metamaterials. 

Maxwell Garnett (MG) approach 

If one of the constituents is regarded as the host material, and the others as 
inclusions, the equation (5.22) can be rewritten as Maxwell-Garnett (MG) equa¬ 
tion [112, 113]: 



(5.23) 


where e m is the dielectric functions of the host material. 

It should be noted that the result depends on whether the first material is em¬ 
bedded into the second or the second material is considered to be embedded in the 
first material. 

Lorentz-Lorenz (LL) approach 

The Lorentz-Lorenz approach is based on the assumption that the host material 
is a vacuum (i.e. = 1) [108, 107]: 



(5.24) 


Effective medium approximation (EMA) or Bruggeman approach 

Assuming that the effective dielectric function acts as a host medium for inclu¬ 
sions, i.e. e = €h, the resulting equation (5.22) can be rewritten with the left hand 
side equated to zero [114]: 



(5.25) 
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Discussion 


One of the approaches is most effective depending on the particular type of 
medium [105]. The MG theory should be applied when constituents clearly may 
be subdivided into the inclusions and the matrix material. The EMA theory works 
best in the presence of molecular mixtures, where a clear subdivision into inclusions 
and the host material is not possible. The LL approach is best suited for porous 
materials. The effective medium theory can be viewed as the first order approxima¬ 
tion, allowing us to obtain the effective permittivity of the host material with the 
inclusions. 

One of the objectives of the presented work is to obtain optical constants for 
various metallic films deposited on the fibre surface with either the chemical or CVD 
technique. Two examples of the him morphology are shown in Figure 5.2. 



Figure 5.2: SEM images of gold him surface morphology, taken at different stages of 
chemical deposition [1], 


The effective optical constant can be obtained in fairly straightforward way using 
the Bruggeman effective medium theory (EMT). An example of EMT application to 
a thin him nucleation and growth on a silica wafer can be found, for example, in [93]. 
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5.3.2 The exact numerical method accounting for the 
interaction between film elements. 

The full numerical simulation of the optical properties of a him can be used as an 
alternative to the simple mixtures models described in the previous section. 

The absorption, and thus the imaginary part of the refractive index, can be ob¬ 
tained if the incident and scattered intensity are known. Next, the real part of the re¬ 
fractive index can be found by implementing the Kramers-Kronig integration method, 
as was described in Section 5.1.1. The idea here was to use the Finite-Difference Time 
Domain (FDTD) method to compute the him absorption properties by illuminating 
it with a short light pulse. The absorbed power can be found by taking the difference 
between the incident and scattered intensities. In order to separate the incident and 
rehected pulses in time domain, either a large simulation area in the space domain 
has to be chosen or the pulse should also be sufficiently short in the time domain. 
However, the pulse should be sufficiently long so that its spectral image is sufficiently 
narrow, allowing the required spectral resolution to be achieved. Therefore, several 
contradictory conditions have to be met. 

We chose “FDTD-method Maxwell solver” from Lumerical Solutions for our sim¬ 
ulations. The software package was designed specifically for electromagnetic held 
simulations and is used extensively. Examples of its applications to the simulation of 
nanopoarticles can be found, for example, in [115, 116]. 

In our research we simulated the rough metal him shown in Figure 5.2. The him 
was approximated by a metallic him substrate in which spherical NP inclusions were 
immersed, as shown in Figure 5.3. The metal him and immersed nanoparticles were 
chosen to have an identical refractive index. Considering the enormous amount of 
time required for FDTD simulation, only a small area with randomly assign rough¬ 
ness was simulated. The sample was assumed to be infinite. This condition was 
satished by surrounding the sample with four planes, transverse to its surface, which 
enforced the periodic boundary condition. The perfect matching layers were placed 
above and beyond the sample, in order to cancel reflection from the computational 
domain boundaries, as shown in Figure 5.3. The pulse was set to incident the him 
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Figure 5.3: Simulation of light scattering by a rough silver him deposited on the glass 
substrate. 

surface from within the glass substrate, which is the case for fibre optical sensors. 
The incident and scattered fields were measured with twelve power sensors (2D rect¬ 
angular plates) detecting the incident held. The hrst six sensors were arranged in 
a 3D cube configuration including the simulation region but excluding the emitting 
antenna. The last six sensors were arranged into the outer cubic shell confining all 
the simulated objects including the antenna. The inner and the outer shells were 
set to detect only ingoing (antenna) and outgoing (scattering) radiation energies, re¬ 
spectively. The absorbed energy was computed by taking the difference between the 
incident and scattered energies. Alternatively, the system geometry can be chosen 
such that the incident and the scattered pulses can be separated in time domain. Re¬ 
peating the measurements for a number of frequencies, in the bandwidth of interest, 
the absorption spectrum was obtained. 

















Chapter 5. Optical properties of materials 


150 


Although the proper 3D him simulation takes considerable amount of time, more 
than 24 h, some useful conclusions can be drawn from a simple 2D simulation. The 
result of light interaction with an array of spherical nanospheres is shown in Fig¬ 
ure 5.4. The obtained simulation results helped to explain the observed anomalies in 
transmission spectra of a metal coated TFBG sensor [2], which was my contribution 
to that work. 



x (nm) 

Figure 5.4: Electric held of S (a) and P (b) polarized light interacting with an array 
of metallic spheres [2], 

It should be further noted that the 2D simulation took less than 1 h, wheres a full 
3D simulation requires more than 24 h. Furthermore, if the goal is to establish polar¬ 
ization or angular dependency, the parametric space should be scanned (by changing 
values of the incident angle and polarization state) which makes this approach pro¬ 
hibitively time expensive. In such a case it would be reasonable to choose a different 
software package capable of running on a supercomputer, or to use a different sim¬ 
ulation technique, for example the discrete dipole approximation, described in the 
following section. 
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5.4 Discrete dipole approximation (DDA) or 

coupled dipole approximation (CDA) method 


The discrete dipole approximation (DDA) is an extremely flexible and powerful 
numerical method used for computing the scattering and absorption of electromag¬ 
netic waves by a target of arbitrary geometry, whose dimensions are comparable or 
larger than the incident wavelength. 

The basic idea of the DDA was introduced in 1964 by DeVoe with application 
to the optical properties of molecular aggregates [117, 118]. At first the retardation 
effects were not included, and treatment was limited to aggregates that were small 
compared with the wavelength. The retardation effects were added in 1973 by Purcell 
and Pennypacker [119]. The accuracy of the DDA method had been confirmed by 
comparing it against the analytical results and experimental measurements [120, 121]. 

The DDA method was popularized by Draine and Flatau. In 1993 they distributed 
discrete dipole approximation open source code DDSCAT [122, 121], which is now 
freely available. DDSCAT is an open-source Fortran-90 software package applying the 
discrete dipole approximation (DDA) to calculate the scattering and absorption of 
electromagnetic waves by targets with arbitrary geometries and a complex refractive 
index. The target may be an isolated object as well as 1-d or 2-d periodic arrays of 
target unit cells. The method can be used to study absorption, scattering, as well 
as electric fields. The description of the latest version DDSCAT 7.2 can be found 
in [123]. At the present day a few other DDA implementations are available. The 
highly optimized computational framework OpenDDA was written in C language 
by J. M. Donald at University of Ireland [124]. Another implementation of DDA, 
called ADDA has been developed through over a period of ten years at University of 
Amsterdam by M. A. Yurkin and A. G. Hoekstra. Their software package is capable 
of running on multiple processors and clusters of computers under Linux Operating 
system [125]. 
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The DDA method implies an approximation of a continuum target (scatterer) by 
an ensemble of polarizable points in a cubic lattice. The polarizable points acquire 
dipole moments in response to the local electric field. Each of the dipoles is driven 
by an incident electric field and by contributions from the other dipoles, hence the 
method is also known as the coupled dipole approximation (CDA) method. These in¬ 
teracting dipoles give rise to a system of linear equations, from which a self-consistent 
solution for the oscillating dipole moments is found, and next the absorption and scat¬ 
tering cross sections are computed. If DDA solutions are obtained for two linearly 
independent states of polarization of the incident wave, then the complete amplitude 
scattering matrix for a given scatterer can be determined as well. An extensive review 
of the DDA method and its applications was given in [122, 126]. The theory behind 
the DDA method was discussed in [127, 128, 125, 129, 130] here we only show the 
idea behind the derivation. 

Let us assume an array of N polarizable point dipoles located at r 3 with each 
one characterized by a polarizability a 3 . The system is excited by a monochromatic 
incident plane wave: 


E loc {r,t)=E 0 j$ ? - ut \ (5.26) 

inducing a dipole moment in each polarizable point: 

Pj = ajEjjoc, (5.27) 

where index j denotes position r) of the jth element and Pj is the dipole moment of 
the )th element. 

Each dipole of the system is subjected to an electric field that can be split in two 
contributions: the incident radiation field E in , plus the field radiated by all the other 
induced dipoles E dp • The sum of both fields is the so-called local field at each dipole: 

Ej,ioc Ejj n T Ej dp • (5.28) 

Now excluding Ejj oc from the consideration by substituting (5.27) into (5.28) we 
obtain 
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— P 

a j 


3 


Ej.dp E n c 


ikri 


(5.29) 


Now let’s consider the field radiated by all the other dipoles Ej^ P - First we use 
the well-known expression for radiation of a dipole located at position fy evaluated 
at point fj 7 ^ f k [40]: 


A kRj 


Ejk 


R 


■jk 


k Tljk X ( 7 I j k X Rk ) + 



ik 


R 


jk 


(3 n jk ■ (n jk ■ Pk ) - Pk) (5.30) 


where Rj k = r} — r kl and nj k is a unit vector pointing from r k point to r} point. 

Taking a close look at equation (5.30) we note that at the given j and k this 
equation simply expresses components of the Ej k vector in terms of linear combination 
of the Pj vector components. Thus we can rewrite this equation in a matrix form: 


Ejk 

where [Aj k \ is a 3 x 3 matrix. 

Now summing all the fields from all k ^ j dipoles we obtain: 


(5.31) 


E j4r = Y} A S*\P (5.32) 

j+k 

Basically we obtained a representation of Ej^ P vector in terms of all the 3 • (N — 1) 
components of all the Pj vectors where j ^ k. 

Finally, (5.29) can be rewritten as: 


—Pj ~ V\A 3k ]Pk = E 0 e 


ikri 


(5.33) 


This is the final equation that has to be solved, consisting of N coupled vectorial 
equations from which N vectors Pj can be found, thus providing a solution to the 
scattering problem. It is more convenient to formulate the problem as a set of scalar 
equations, therefore we can simply rewrite the above equations as: 


[A}P = E m 


(5.34) 
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For N dipole elements the [A] is a 3 N x 3N symmetric matrix, P = (Pi , Pn) 
and E in = (P 0 e*^ ri , E a e l ^ rN ) are 3N vectors. 

Once the polarizations Pj are found, the extinction, absorption and scattering 
cross section can be expressed as follows [131]: 

farb JL 

°ext = TTT-^ Im i E hn ' Pj], ( 5 ‘ 35 ) 

f-Jin 

J 

the extinction cross section is computed from the forward-scattering amplitude using 
the optical theorem [131]; 


' j j j AC i j 

Cabs = rrr-^ Im [ E j,in ■ aE i n l ( 5 - 36 ) 

| -C^in I 

the absorption cross section is obtain by summing rate of energy dissipation by each 
of the dipoles; 


r< 

^ am. 


ink 
I E- I 2 

| in | 


« n 

J £[3 


D —ik(nrj) 


dVt, 


(5.37) 


where n(Q) is a unit vector in direction of scattering and dil is the solid angle element. 
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5.5 Conclusion 

In this chapter we have reviewed the optical properties of materials, briefly dis¬ 
cussed experimental methods which are used to determine the optical constants n and 
k. We have chosen parameters for eleven metals commonly used in optics and plotted 
their dispersion curves, which are extensively used in the following simulations. 

The morphology of various films can be modeled with the simplest mixtures model. 
We are going to use the Maxwell-Garnett model to compute the effective refractive 
index of rough films, and the Lorentz-Lorenz model to determine the refractive indices 
of liquid mixtures. 

The properties of nanoparticles are discussed in the next Chapter, where we will 
obtain absorption coefficients and deduce the effective optical constants of nanopar¬ 
ticles. The obtained effective permittivity of nanoparticles can be considered as the 
permittivity of inclusions in the Maxwell-Garnett model. Thus it should be possi¬ 
ble to characterize metallic films of various roughnesses as well as films consisting of 
nanoparticles. 

A more accurate result can be obtained if the radiative dipole-dipole interactions 
between the film elements are considered. Unfortunately the FDTD method is pro¬ 
hibitively time expensive, although we did reveal some useful results. 

Alternatively we might consider the DDA method. For the sake of simplicity, 
properties of a particular particle can be found with other methods, such as the 
Mie theory discussed in the following Chapter, and then used to define the starting 
points of the DDA method. Therefore all the elements of nanoparticle-based coating 
become intrinsically interconnected, thus allowing us to account for a collective effect. 
However, for sparse coatings the analysis can be solely based on optical properties of 
one single nanoparticle. 




Chapter 6 


Optical properties of nanoparticles 


In this chapter we discuss optical properties of metal nanoparticles. Next we 
present results of our simulations. The results of simulation would allow us to deduce 
parameters for an optimal nanoparticle based coating enchaining the TFBG sensor 
sensitivity. 

6.1 Review 

6.1.1 Characterization of light scattering by particles 

We start by introducing parameters which are usually used for particles character¬ 
ization. Let us consider a single particle illuminated by a plane wave with intensity 
Ii nc , than the total power scattered by this particle is W sca , defined as: 

ll.scn Csea line •: (6.1) 

here C sca is the scattering cross section, with dimensions of area. 

Particles can also absorb electromagnetic radiation. The rate of absorption W a f )S 
should also be proportional to the incident intensity: 

^'V a bs = CaJasIinci ( 6 . 2 ) 

where C a b s is the absorption cross section. 
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The sum of the scattering and absorption sections is called the extinction cross 
section: 

C e xt = Csca + Cabs■ (6-3) 

In practice, the cross sections are usually normalized by the particles area pro¬ 
jected onto a plane perpendicular to the incident beam, and defined as Q ex t, Q SC a, 


polarization-dependent properties 


Particles can be thought as miniature polarizers and retarders. Assuming that 
the incident field is polarized, the polarization of the scattered field would depend on 
the nanoparticle shape, size and the direction of scattering. Even unpolarized light 
can become partially polarized upon particular direction of scattering. 

The polarization properties of a particle can be described, for example, with Jones 
matrix. Assuming that the incident wave, with the wavenumber k, propagates along 
the 0 axis, the general relation between incident and scattered fields can be written 
in the following from [132]: 

E sca = e —^[S}E mc . (6.4) 

For the known direction of incident light and a particular chosen direction of the 
scattered light, the scattering plane can be introduced, spanned by these two direc¬ 
tions. Now fields can be decomposed into two orthogonal components, one parallel, 
the other perpendicular to the scattering plane: 



0 ik(f—z) 

—ikr 


Su s 12 
S21 S22 



(6.5) 


The elements Sjk of the scattering matrix [S] (or Jones matrix) are complex-valued 
functions of the scattering direction defined, for example, by 9 and (ft in polar coordi¬ 
nate system. 
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The scattering Mueller matrix can also be introduced, connecting the Stokes pa¬ 
rameters, which can be measured directly [132]: 


(l \ 
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here 


Q = ^(e,e;-e iE1) , 


(6.7) 


U = 2 k ~^ E - + 

V = 2^ (E « E1 ~ E-lE;)- ( 6 - 8 ) 

It can be noted that all the parameters I,Q,U,V are real numbers and can be mea¬ 
sured by optical means with quadric detector, without the necessity to measure the 
phase of the signal, unlike in the case of characterization by the Jones matrix. This 
subject was discussed in more detail in Chapter 4. 

An example calculation can be found in reference [133], where polarization prop¬ 
erties of light scattered by four different homogeneous particles: a prolate spheroid, 
an oblate spheroid, a finite cylinder and a bisphere with touching components were 
considered. The exact analytical solution exists for such scattering problems, which 
made it possible to verify T-matrix and DDA methods techniques [133]. The particle 
was illuminated by a plane harmonic wave propagating along the positive z-axis. The 
scattering matrix [5] = Sjj was computed as a function of the scattering angle. It 
was shown for all scattering angles the elements S 13 , Su, S 23 , S 24 , S 31 , S 32 , S 41 , S 42 
vanish and Sn = S 22 , *S'i 2 = S 21 , S 34 = — S 43 , S 33 = S 44 . Therefore only the remain¬ 
ing elements Sn, S 21 , S 33 , S 43 were considered. A strong polarization dependence 
was revealed, especially in the case of prolate spheroid. Another example, revealing 
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a strong polarization dependence of optical absorption can be found in [134], where 
a prolate spheroid with an aspect ratio of 1.6 and a major axis of 8 nm embedded in 
silica was illuminated with linear polarized plane wave. 

The extensive study of polarization properties of light scattered by spherical par¬ 
ticles can be found in the original Mie’s paper [135] as well as in [136, 137], where it 
was shown that for a particle smaller than lOOnm the scattered light behaves accord¬ 
ing to the Rayleigh scattering law and has the maximum polarization at scattering 
angle 6 = 90°. Even if particles are illuminated by unpolarized light, the scattered 
light becomes partly polarized. The scattering patterns of a spherical particles were 
also considered in [136], where computation method was based on numerical evalu¬ 
ation of Mie’s coefficients. It was shown that the maximum of polarization moves 
to 6 = 120° with the increase in the size of a particle. For particles with sizes larger 
than lOOnm the degree of polarization diminishes rapidly. Light scattered sideways 
is always linearly polarized, regardless of the particle size [136, 137]. 

We conclude, that even spherically symmetric particles have nontrivial polariza¬ 
tion properties which depend on the scattering angle, particle size and material prop¬ 
erties from which the particle is made. Even a small variation in material absorption 
is causing significant change in the polarization properties of scattered light. 

Size dependent effects 

The absorption and scattering phenomena can be viewed as a size-dependent pro¬ 
cess. The first extensive review of this subject was given in Mie’s original paper [135], 
where absorption and scattering spectra were calculated for gold nanoparticles with 
various sizes. 

The classical Mie’s approach can be applied with no limitations to nanoparticles 
larger than 20 nm [138]. For smaller particle the material properties would differ 
significantly from measured bulk material properties. The imaginary part of the 
dielectric function, which is responsible for the absorption, increases by a factor of 10 
compared to the bulk value for particles of 2.5 nm whereas for 20 nm particles the 
difference is only of a factor 1.5 — 2, depending on wavelength [138]. 
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For a particle of a few nanometers in size the mean free path of conduction elec¬ 
trons is limited by the particle boundary, which interferes with the electron relax¬ 
ation time and leads to the increase in absorption at a narrow wavelength range. 
This effect can be taken into account by introducing an extra damping term into the 
Drude-Sommerfeld free-electron model. 

The experimental investigation performed by Soonnichsen [138] confirmed that 
the quantum treatment is not necessary for metallic particles larger than 20nm. It 
was also shown that the collective oscillation dephasing is adequately described by 
single electron dephasing, following from the classical solution. The surface scattering 
effect and chemical damping by the outside medium can be negligible as well. Thus 
the bulk material properties can be successfully applied to simulation of particles with 
sizes larger then 20 nrri. 

Shape and substrate dependent effects 

As was mentioned previously, the exact analytic solutions to the scattering prob¬ 
lems is known only for a few simple geometries. Numerical methods, such as discrete 
dipole approximation, have to be applied in the general case. 

The size dependent scattering and absorption properties of non-spherical particles 
were studied, for example, in [139]. Discrete dipole approximation (DDA) method 
was chosen for simulations of absorption and scattering spectra of gold spheres, cubes 
and rods, ranging in size from 10 to 100 mu . The problem of light scattering by silver 
particles with various geometries (spheres, cubes, truncated cubes) was also reviewed 
in [134], First the polarizability was calculated using the Clausius-Mossotti relation 
with dielectric function of bulk silver. It was observed that the optical response below 
325 nm is independent of a particle shape and is defined by intra-band transitions in 
silver. 

The optical properties of particles are also strongly dependent on the substrate 
properties. For example, if a nanoparticle is deposited on a substrate, the substrate 
can be considered as an array of dipoles interacting with the nanoparticle [140]. 
Alternatively, the model can be simplified by replacing the substrate with an induced 
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image of the nanoparticle, thus dipolar interaction between the particle and its image 
can be considered instead [134]. However such dipolar interaction leads to the strongly 
inhomogeneous held, especially when the particle is close to the substrate. In such a 
case multipolar interactions have to be considered [141]. 
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6.2 Simulation of light scattering by small 
particles 

The problem of light scattering by a spherical particle made of an arbitrary mate¬ 
rial, including absorbing materials such as metals, can be solved exactly by analytical 
methods. The exact solutions are also known for spheroids, or infinite cylinders. 

The solution is usually obtained by expanding the incident, scattered, and internal 
fields in a series of vector spherical harmonics. The expansion coefficients are found 
from the boundary condition ensuring continuity of tangential components of the 
electric and magnetic fields across the surface of a particle. Observable quantities are 
expressed in terms of the coefficients. 

An interesting overview of the early work was given in [142, 137, 143], an ex¬ 
tensive Mie’s theory review can be found in [136, 143]. As early as in 1863 Alfred 
Clebsch (1833-1872) published a general solution of the elastic wave equation in 
terms of the vector wave functions [144] and Ludvig Lorenz (1829-1891) completely 
solved the problem in terms of the ether theory. Lord Rayleigh (John Strutt) (1842- 
1904) introduced his theory of elastic light scattering by small dielectric particles in 
1871 [145, 146]. The particles were considered to be much smaller than the wave¬ 
length of the light so that the general electromagnetic problem can be reduced to 
electrostatic problem of an isotropic, homogeneous, dielectric sphere in a uniform 
held. 

The interaction of electromagnetic waves with a perfectly conducting sphere was 
reviewed by Thomson [147]. The solution was generalized by Hasenorl [148] who 
introduced a non-zero conductance for metals. His paper actually already gives the 
full solution contained in Mie’s later publication [137]. Ehrenhaft [149] gave a rigorous 
treatment of the scattering of light by small absorbing spheres, which is an even more 
elegant than Mie’s work [137]. A similar treatise was given by Debye [150], who 
provided an alternative approach by utilizing two scalar potential functions. 

Finally in 1908 Gustav Mie presented his solution to the problem in recognizably 
modern notation, and presented a comprehensive comparison of experimental results 
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with theoretical hirelings, providing many numerical examples [135]. Although in his 
calculation Mie was using only three terms in the expansion series it was sufficient 
to completely explain all the optical phenomena observed until then. All Mie’s con¬ 
clusions published in his paper from 1908 remain valid. Today the problem of light 
scattering by a homogeneous isotropic sphere is firmly attached to the name of Gustav 
Mie, although Mie was rather the last one who solved this problem. 

Mie had to limit his results to particles smaller than 200 nm due to the limits 
imposed by calculation. The calculations were performed by hand, and only three 
first terms in the infinite series were considered. He also limits his finding to spherical 
particles, although he knew that ellipsoidal particles can be treated in the similar 
way [137]. Today, with availability of computers, these limitations can be easily 
removed. 

The exact solutions to Maxwell’s equations are known only for special geometries 
such as spheres, spheroids, or infinite cylinders, so approximate methods are required 
for the general case. An excellent review of various methods can be found in [151, 
152, 153, 154, 129, 127]. Mie’s theory is usually used as reference, to validate other 
methods. 

6.2.1 Analytical solution to the problem of electromagnetic 
wave scattering on spherical particles 

In this section we briefly outline the steps taken to solve the problem of elec¬ 
tromagnetic wave scattering on small spherical particles. An excellent review of the 
subject can be found in [126, 127, 128, 132], 

We start with formulation of the problem. Let us assume that an incident elec¬ 
tromagnetic field E inc , H inc interacts with a particle that occupies a bounded region. 
The total fields E tot , H tot in the surrounding medium are equal to the sum of the 
incident and the scattered fields: 
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We are interested in a free space solution outside the particle boundaries, therefore 
we assume that there are no sources in space. Hence, Maxwell’s equations can be 
reduced to the vectorial Helmholtz equations, as was shown in Chapter 2: 


(V 2 + k 2 {r))E{f) = 0, 
(V 2 + fc 2 (f))#(f) = 0, 


( 6 . 10 ) 


where 



( 6 . 11 ) 


Now let us consider the boundary condition: (i) the tangential components of 
the total fields have to be continuous across the particle surface; (ii) the radiation 
condition, requiring that the tangential components of the electric and magnetic fields 
have to approach zero at rate 1/r as the distance r from the origin approaches infinity. 

The solution to the vector Helmholtz equations (6.10) can be searched in terms 
of scalar basis functions, obtained as a solution to the scalar Helmholtz equation: 


(V 2 + k 2 (‘r))ip('r) = 0. 


( 6 . 12 ) 


In spherical coordinates equation (6.12) becomes: 



1 d 2 , 1 8 ( . 

r dr 2 ^ + r 2 sin(0) 89 


( sia{e) T) 



Separating the variables with ansatz 


^M,0) = i*(r)0(0)$(0), 


(6.14) 


the three ordinary differential equations are obtained: 


The solution to the radial equation in (6.15) can be expressed in terms of special 
functions 



R(r) = uf\ 


( 6 . 16 ) 
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where 

u\ 1 ^ = Ji(r ) are spherical Bessel functions, 

/o\ 

u\ ; = Ni(r) are spherical Neumann functions, 

= H[ l> (r) = Ji(r) + iNi(r ) are spherical Hankel functions of the first kind and 
■u[ 4) = Hj 2> (r) = Ji(r) — iNi(r) are spherical Hankel functions of the second kind. 
The solutions to the polar equation in (6.15) are the associated Legendre functions: 

0(0) = pr(e), (6.17) 


where m = 0, ±1,.... ± l. 

Finally the solutions to the azimuthal equation are simply harmonic functions: 

<f>(0) = e irnrt> . (6.18) 

Thus the solution to the scalar Helmholtz equation in spherical coordinates can 
be written in the following form: 

(6.19) 

Using the above scalar basis function we can create vectorial basis functions: 

££(5) = v ■«>,£• 5), 

«£(3) = (vxM« (6.20) 

Here a is some constant vector or the position vector a — f, in such a case func¬ 
tions (6.20) satisfy the vector Helmholtz equation in spherical coordinates (6.10) [127, 
132], 

However, we assumed that Maxwell’s equations are divergence-free , i.e. the 
particle is immersed in a charge free medium, but the functions L^(a) are not 
divergence-free and therefore should be excluded from consideration. Hence, the 

incident, scattered and internal electric fields can now be represented in terms of 
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spherical vector wave functions: 



L l 

- EE 

Z—0 m=—l 




L l 

= EE 

1=0 m=—l 



Eint^T ) 

L l 

= EE 

1=0 m=—l 

( c ZnM^(kf) + cf m A^(A;r)) . 

(6.21) 


Here E int defines the field inside the particle and k is the wavenumber inside the 
particle. 

To satisfy the boundary condition at the origin a non-singular at the the origin 
spherical vector wave basis functions of the first kind have to be used. The scattered 
field has to satisfy the radiation condition, therefore it can be expanded in spherical 
vector wave functions of the third kind. 

Considering the continuity of the tangential field components across the surface 
of the particle, the boundary condition can be set: 

x (Ei nc (k 0 r ) T E sca {k 0 r ) Ei n f (kv ) ) 0, 

n+ x ( H inc (k 0 r) + H sca (k a r) - H int (kr)) = 0. (6.22) 

Here n + denotes the outward pointing vector, normal to the boundary surface of the 
particle. 

The boundary conditions yield a system of linear equations, from which the un¬ 
known expansion coefficients f3 ^ n , f3[* m and , c^ m can be expressed in terms of the 
known expansion coefficients a^ m , characterizing the known incident field. 

For spherical particles this system of linear equations can be inverted analytically, 
and leads to the well-known Mie’s solution [135, 150]. 

Observable quantities such as scattering C sca , absorption C a b s and extinction 
C ex t = C sca + Cabs cross sections can be expressed as follows [132]: 
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C ex t — 


2ir °° 

— ^^(2 1 + 1 )Re(ai + bi ), 

i=i 
2vr 


^ = -u + z +u*)■ 

1=1 


(6.23) 


Here 


ai = 


bi = 


ml)i(nx)'il)' l {x) — ijji(x)'ijj , l (nx) 
r mjji(nx)^ , l (x ) — ^i{x)il}[{nx) 


(6.24) 

i>i(nx)£[(x) - n£i(x)^[{nx) ' 
with size parameter x = ka, where a is the radius of the sphere and k is the wavenum¬ 
ber of the incident light in the surrounding medium, and £; are the Riccati-Bessel 
functions, and n = yfe is particle refractive index. 

Alternatively the Mie solution can be derived with help of Debye potentials [110]. 

The described procedure can be applied in any other coordinate system in which 
the scalar Helmholtz equation becomes separable. Hence various particles with non- 
spherical symmetry can be considered as well [155, 156]. However, the basis functions 
might no longer be orthogonal and a resulting system of linear equations can be only 
inverted numerically. 


6.2.2 The quasi-static approximation 

In this section we review the quasi-static approximation, which can significantly 
simply analysis of non spherical particles in the following section. 

For a particle with size significantly smaller than the wavelength of the incident 
light, the electric held of the light can be viewed as an uniform and static. Thus the 
interaction between a small particle and light can be described in terms of electrostat¬ 
ics rather than electrodynamics. The electron cloud of the particle can be viewed as 
displaced by the electric held, hence a restoring force arises from Coulomb attraction 
between electrons and nuclei, resulting in electron cloud oscillation. 

If a homogeneous, isotropic sphere is placed in medium with a different permittiv¬ 
ity, subjected to a uniform static electric held, a charge will be induced on the surface 
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of the sphere. The electric fields inside and outside the sphere (of radius R) can 
be described by Laplace’s equation in terms of scalar potentials 0i(r, #) and 0 2 (r, 9) 
written in spherical coordinates [132, 40]: 


V 2 0i(r < i?,d) = 0, 

V 2 0 2 (r > R,9) = 0, (6.25) 

At the boundary between sphere and medium the potentials must satisfy conti¬ 
nuity condition: 


MW) 

—<9 r 0i (R,B) 


02 (R,0), 
—d r 4> 2 (R,9). 


(6.26) 


Thus we have the partial differential equation and two boundary conditions. In 
addition, by requiring the held at infinity to be equal to the external homogeneous 
held E 0 , the solution can be written in the following form: 


0i(r,0) =- E 0 rcos(9), 

ei + ze 2 

Mr, 9) = ~E 0 r cos (9) + R 3 E Q £l ~ 62 C ° S f } . (6.27) 

ei + 2e 2 r z 

Now superposing the above two potentials and comparing the result with the 
potential of an ideal dipole: 


0M) = 


1 p ■ r 1 p cos 9 


47re 2 r 3 47re 2 r 2 
The solution (6.27) can rewritten in the following form: 


(6.28) 


V = ^2 aE 0 , 


a = 47ri? 3 


€l - 62 

ei + 2e 2 ’ 


(6.29) 
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where a(R, ei, e 2 ) is the polarizability, which depends on the particle size and material 
permittivity, e\ is the material permittivity, from which the particle is made and 62 
is the medium permittivity. 

The above result states that a sphere in an electrostatic field is equivalent to an 
ideal dipole, with momentum p = 62 aE 0 . Therefore the problem of electromagnetic 
wave scattering by a sphere can be reduced to a well studied problem of dipole 
radiation. The cross section for absorption and scattering of field by an ideal dipole 
exposed to a plane wave source is given by the following relations [132, 40]: 


C^s = k$s(a), 

C S ca = |h«| 2 , (6.30) 

bn 

here k — n 2 ^ is the wavenumber of incident wave in the medium outside the particle. 
However it is more convenient to use in practice unitless parameters called efficiency: 


Qabs 


Q 


sca 


c, 


abs 


nR 2 ’ 


nR 2 ’ 


(6.31) 


6.2.3 The ellipsoidal shape particles 


The analytical solution for ellipsoidal shape particle can be found in a similar 
way as in the case of spherical particle by introducing ellipsoidal coordinates and 
exploiting particle symmetries to solve Laplace’s equation. The derivation can be 
found for example in [132, 40]. Here we review the quasi-static approximation instead. 

As before we solve the problem of electric charge distribution on a particle of a 
given shape subjected to a uniform static electric field. The result can be written 
in the similar form as in the previous section. For a cigar-shaped ellipsoid with two 
equal length principle axes, defined by b and a > b (such spheroid is called prolate) 
we have: 


p = e 2 aE 0 , 

_ 4vr 2 ei - e 2 

CL X — “ 7 ~ CIO --— -—, 

3 e 2 + L x (e\ — e 2 ) 


(6.32) 
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The material permittivities for the particle and the surrounding medium are defined 
as 6i and e 2 , respectively. The parameter L x is called depolarization factor and 
depends on the relative orientation of the long axis a of prolate spheroid and the 
polarization direction of incident light, which is incident transversely to the plane 
on which particles are deposited. For incident light with the polarization direction 
parallel to the long axis of the prolate ellipsoid x = a the depolarization factor is 
defined as: 


L n = 


2e 


ln 


1 + e 


- 1 


(6.33) 


and for the polarization direction perpendicular to the long axis x = b the depolar¬ 
ization factor Lb can be defined in terms of L a [132]: 


L b 


1 

2 


(1 - L a ) , 


(6.34) 


here parameter e = %Jl — r\ 2 is the eccentricity, with the aspect ratio = K 

The cross sections C a b s , C sca , C ext defined in the previous section by equation (6.30) 
can be applied without modification, but the efficiencies in (6.31) should be modified 
by replacing the particle radius with the product of its principle axes: 


Qj 


7r (a6 2 )i 


(6.35) 


6.3 Conclusion 

At this point we know how to define optical properties of metals, how to simulate 
properties of separate nanoparticles, either exactly or with help of the quasi-static 
approximation, and how to model films with different morphology. In the next Chap¬ 
ter we apply this knowledge to design the optimal coating for TFBG based refractive 
index sensors. 







Chapter 7 


The optimal parameters for a 
nanoparticle-based coating 


In this chapter we study optical properties of metal nanoparticles embedded in a 
homogeneous ambient medium, with particular attention to the absorption properties 
of nanoparticles. The absorption properties allow us to compute the effective refrac¬ 
tive index of a medium by applying the Kramers-Kronig relation and connecting the 
imaginary and real part of the refractive index. Next we discuss the choice of an 
optimal nanoparticle-based coating. 

7.1 The idea behind sensitivity enhancement 

Although the optical properties of bulk materials are known, the optical properties 
of a particle can not be directly deduced. The particle size and its shape introduce 
additional resonances and hence additional absorption. The absorption peaks are 
not only influenced by the shape and size of the particle but also by the external 
mediums refractive index. Hence, if the refractive index of the external medium is 
changed the particles resonances will change their locations as well. Thus a coating 
layer can be created with optical properties defined by the refractive index of an 
external medium. In such a material not only the imaginary part of the refractive 
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Figure 7.1: The shift in absorption peak (red) and the corresponding shift in real 
refractive index of a medium consisting from nanoparticles. 

index becomes sensitive to the environmental changes, but also the real part through 
the Kramers-Kronig connection, as shown in Figure 7.1. 

We are mainly interested in the real part of a materials effective refractive index, as 
it defines the positions of resonances, which can be observed in the TFBG spectrum, 
and hence the change in external medium can be easily detected. 

The connection between high quality factor resonances of TFBG structure and 
low quality factor resonance of nanoparticles is schematically shown in Figure 7.2. 

7.2 Simulation of optical properties of metal 
nanoparticles 

7.2.1 The quasi-static approximation 

Let us first simulate optical properties of small spherical nanoparticles made of 
various metals. 

Following the results from Chapter 5 we obtain complete characterization of opti- 
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Figure 7.2: The schematic representation of connection between high Q resonances 
of TFBG and low Q resonances of nanoparticles. 


cal properties of 11 metals in the wide range of frequencies (from uj — 0 up to 10 — 20 
eV ). The metals were described phenomenologically by Lorentz-Drude model, in 
which parameters of six oscillators were fitted for the best consistency with experi¬ 
ments [5]. 

We start with the quasi-static approximation, which is valid for the problem of 
light scattering on spherical particles with the particle size significantly smaller then 
the wavelength. The wavelength of interest is A = 1.5 im , hence the the approxima¬ 
tion should yield the correct result for particles about 100 nm in diameter. 

Considering the Lorentz-Drude model from Chapter 5 and using equations (6.29)- 
(6.31) the absorption and scattering efficiencies can be plotted for various metals as 
shown in Figure 7.3. 

The next Figure 7.4 shows the absorption and scattering efficiencies of gold nanopar- 
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tide immersed in solutions with various refractive indexes. 



Frequency ro (eV) 

Figure 7.3: Simulated absorption efficiency of 30 nm spherical nanoparticle made of 
Au, Ag, Cu and A1 metals, as a function of photon energy (the graphs for Ag and A1 
were divided and multiplied by 10, respectively). 




Figure 7.4: Simulated absorption and scattering efficiency of 30 nm gold spherical 
NP, immersed in media with various refractive indexes. 


From the above results we can draw the following conclusions: 
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1. From equations (6.29) and (6.31) we note that the absorption Q a bs{u) and scat¬ 
tering Qsca(u) curves are independent of the particle size, hence the position of 
plasmon resonance is determined only by the dielectric functions of the particle 
6i and the surrounding media e 2 . The position of the resonance is given by the 
pole of the polarizability in equation (6.29), i.e. when + 2e 2 is approaches 
to zero. This pole is commonly referred as a polarization mode or plasmon 
resonance. 

2. The position of the absorption resonance is different for different metals. As 
can be seen from Figure 7.3 for Au, Ag and Cu metals the absorption resonance 
is located in the visible band, wheres for A1 it is located at IR band, not far 
away from the interband absorption peak. We can also conclude that the A1 
coating might look promising with regards to sensing application in IR band. 

3. Although we used a first order approximation, the quasi-static approximation, 
valid only for a particle with size significantly smaller than the wavelength, we 
came to the correct qualitative as results reported in the literature. The ob¬ 
tained absorption spectra shown in Figure 7.4 are in a good agreement with ex¬ 
perimental measurements, predicting the red colors of light transmitted through 
a suspension of fine gold particles. This effect happens due to the strong inter¬ 
band absorption of the shorter wavelength by the gold. 

4. With the increase of a particle size the electronic cloud displacement can ex¬ 
hibit high-multipolar charge oscillations, thus several poles in the polarizability 
function can appear and significantly affect the absorption spectrum. The quasi¬ 
static approximation can no longer be applied, hence the exact Mie’s theory has 
to be used, as shown in the next section. 
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7.2.2 The exact solution 

As was discussed in Section 6.2.1 the problem of light scattering by a spherical 
particle can be solved exactly. The absorption, scattering and extinction efficiencies 
can be written in the form of infinite series [132]: 


Qext(x) = 7^^(2I + l)%(a:)+6i(a:)), 

' ' i=i 

2tt N 

= ^£( 2 ' + i)(M*)l J + IM*)l 2 ). 

Qabs([•&') R'ext ( x ) Csca(*£)j 
with ai(x) and bi{x) given by: 


(7.1) 


ai(x) = 
bi(x) = 


rmjji(mx)d x ^i(x) - ipi(x)d x ipi(mx) 
rmj}i(mx)d x ii{x) - £i(x)d x il>i(mx) ’ 
i/ji(mx)d x ^i(x) - rmpi(x)d x 'ipi(mx) 


(7.2) 


^i(mx)d x ^i(x) - m£i(x)d x i/ji(mx) ' 

Here x = k a R is the size parameter, dependent on the radius R of the particle, 
the wavenumber k Q of the incident light ( k 0 — ^ ) and the relative complex refractive 
index of the sphere m = n partic i e /n rner i lur n■ The wavenumber k corresponds to the wave 
in the surrounding medium. 

The Riccati-Bessel functions ipi{x) and £i(x) are defined as follows: 


ijji(x) = 

Zi( x ) = 



(7.3) 


where Hl(x) = Ji(x) + iYi(x) is the Hankel function of the first kind, and Ji(x) and 
Y[(x) are Bessel functions of the first and the second kind. 

The number N in equation (7.1) is the number of terms in Mie series. In order 
to apply these equation in practice the series must be limited to a finite number of 
terms N max . Analysis of the convergence behavior of these series reveals that the 
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sufficient number of terms, dependent only on the size parameter x, is determined by 
the following equation [157]: 

Nmax = X + 4:X^ + 1 (7.4) 

The relation (7.4) applies to all refractive indices since the series (7.1) convergence 
is determined entirely by Bessel functions of x alone, as can be seen from (7.2). Mie 
himself considered only three terms as all calculations were done by hand. At the 
present moment it’s possible to compute hundreds of terms in several minutes on a 
typical PC. 

The terms in the Mie series were computed with the help of Mathematica software 
from Wolfram Research. The code is presented in Appendix 9, and results for the 
absorption, scattering and extinction efficiencies are shown in Figure 7.5. 



Figure 7.5: The absorption, scattering and extinction efficiencies as functions of the 
size parameter x , for a particle with the relative refractive index m = 5 + j0.4. 

Next, considering optical properties of metals we can compare absorption efficiency 
for 30 nm spherical nanoparticle made of various metals, as shown in Figure 7.6. The 
figure also shows the results based on the quasi-static approximation. 

It can be noted that the proximate and the exact results are in good mutual 
qualitative correspondence, although quantitative differences are also noticeable, es¬ 
pecially for nanoparticles made from silver. This confirms the well known fact, which 
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Frequency 00 (eV) Frequency © (eV) 

Figure 7.6: The absorption efficiency of 30 11111 spherical nanoparticle made of various 
metals ( Au, Ag, Cu and A1) as a function of photon energy. The result were obtained 
with use of the exact Mie’s theory and with the quasi-static approximation theory. 

states that the optical properties of small metal particles are predominately defined 
by the bulk material properties from which the nanoparticles are made. 

Next we will consider the influence of the refractive index of the surrounding 
medium on the absorption and scattering spectra of gold nanoparticles. The results 
are shown in Figure 7.7. Comparing these results with the quasi-static approximation 
Figure 7.4 we can conclude that the exact Mie solution predicts a significantly different 
response, with a smaller correlation between the change in the external refractive 
index and the absorption and scattering spectra. 

The most striking difference is observed when the absorption and scattering effi¬ 
ciencies are plotted against the particle radius, as shown in Figure 7.8. In the case of 
the exact solution a number of local resonances is observed, where in the case of the 
quasi-static approximation only a single resonance is predicted. The single resonance 
is defined by the approximation equation (6.30) where Q abs R and Q abs R 4 , in 
accordance with the Rayleigh scattering theory [145, 146]. 

Our goal here is to find the optimal set of system parameters, such that the local 
plasmon resonances of nanoparticles can be effectively coupled to the cladding modes 
















Chapter 7 . The optimal parameters for a nanoparticle-based coating 179 



Frequency © (eV) 



Figure 7.7: Absorption and scattering efficiency of 30 nm gold spherical nanoparticlc, 
immersed in media with various refractive indexes. 



Figure 7.8: The size dependence of the scattering and absorption efficiencies of gold 
NP illuminated by electromagnetic wave at A = 560 nm. The particle is assumed to 
be in the medium with refractive index n — 1. 


excited in the TFBG sensor. We have a particular interest in the absorption efficiency 
Qabs of nanoparticles as it affects the effective refractive index of the surrounding 
medium detected by the sensor. With this in mind we can plot Q a b s as a function of 
various parameters, and thus study the parametric space. 
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First, let us study the dependence on the particle radius and the material optical 
properties from which particle is made. Figure 7.9 shows the absorption efficiency 
Q a b s as a function of particle radius and incident photon energy. The distinct location 
of the plasmon resonances is intrinsically connected with the material from which the 
particle is made. It is also should be noted that as the particle size increases higher 
modes of plasmon oscillation can occur, such as quadrupole modes and higher order 
modes. However such higher order modes are not always excited, for example if the 
particle is illuminated by a plane wave, therefore the energy absorption might not be 
observed in the spectrum. Such modes are often referred as dark plasmon modes. 

We can summarize our findings as follows: 

1. Both theories: the approximation quasi-static theory and the exact Mie theory 
are predict similar qualitative results, stressing the importance of the metal 
dielectric function from which the particle is made. However the predictions 
vary significantly for large particles. 

2. From Figure 7.9 it can be seen that there is an optimal particle size for each given 
wavelength, at which the absorption peak has the greatest value. This result 
was rather unexpected, but a similar conclusion can be drawn from Figure 7.8, 
obtained in several other papers as well. 

3. The material choice is extremely important. If we seek to excite plasmons in 
the visible band, a nanoparticle made of gold is a good option. However, as 
is seen from Figure 7.9, copper can also be a good choice for a relatively large 
nanoparticle. If the objective is to obtain a resonance in the IR band aluminum 
might be a reasonable choice, although the resonance is extremity weak. 



Particle Radiun (nm) Particle Radiun (nm) 
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Figure 7.9: The absorption efficiency of gold, silver, aluminum and copper particles as 
a function of particle size and incident photon energy (plotted in logarithmic scale). 
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7.2.3 The shape effect. Ellipsoidal nanoparticles. 


In the previous section we have established that the quasi-static approximation 
theory can provide satisfactory results in the case of small particles. This conclusion 
was drawn by comparing the exact scattering theory with the quasi-static approxima¬ 
tion. Therefore we can use the quasi-static approximation to investigate the influence 
of a particle geometrical shape on its optical properties. 

We have a particular interest in randomly orientated silver nanorod particles, 
shown in Figure 7.10, as some interesting experimental result were reported [3]. 
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Figure 7.10: SEM image of the fibre surface coated with gold nanorod particles 


Such nanorod particles can be treated as prolate spheroids, with a significant ratio 
between its principle axes a and b. Random orientation can be considered simply by 
taking the average polarizability [158]: 

1 2 . . 

ol = —OL a + -a b , (7.5) 


here a a and a b are polarizabilities along a and b principle axes, respectively. The 
cross sections for absorption and scattering processes can be found with help of equa¬ 
tions 6.30 and 6.31, as was discussed previously. 
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The simulation result is shown in Figure 7.11. It can be observed that silver and 
gold particles have distinct resonances in the visible band and look promising for the 
further investigation. 




Figure 7.11: Comparison of the optical absorption and extinction efficiencies of prolate 
ellipsoids made from various metals with the principal axes b = 30, a = 100 nm. 

We next show that by changing the aspect ratios between principal axes of an 
ellipsoid particle we can tune the absorption resonance into the operational range of 
TFBG sensor. The figure 7.11 shows the absorption efficiency of randomly orientated 
prolate silver particles with various aspect ratios. 

The properties of the TFBG sensor coated with elongated nanoparticles might be 
simulated in the usual manner, except that we have to account for the film anisotropy, 
particularly for the difference in optical properties of the film seen by the tangential 
and transverse polarized electric field. The difference arises from the fact that the 
position of the absorption resonances of an elongated nanoparticle is a polarization- 
dependent function, hence the real part of the effective dielectric permittivity function 
is also polarization-dependent. It is convenient to introduce effective dielectric per¬ 
mittivities along the tangential and transverse axes associated with the film planes, 
as shown in Figure 7.13. 
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Figure 7.12: Comparison of the absorption efficiency for prolate silver nanoparticles 
with different aspect ratios between the principal axes r/ = - 



Figure 7.13: The asymmetry of the coating made of elongated nanoparticles, a) The 
transversely polarized electric held E p encounters e p dielectric permittivity and b) the 
tangentially polarized electric held E$ sees the e $ dielectric permittivity. 


The basis functions of a cylindrical waveguide coated with elongated nanoparticles 
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can be found from the modified equation ( 2 . 21 ): 


^ + ~ p dp + d n< 7>) 



+ e p kl + { lnep )"-/? 2 



+ e <A'o — P 2 _ 



0 , (7.6) 


here e p and are the dielectric permittivities as seen by the transversely and tan¬ 


gentially polarized electric fields E p and E^, respectively. 


Our main interest is to increase the TFBG sensor sensitivity operating the in the 
IR band. We can expect that by coating the sensor with silver or gold prolate parti¬ 
cles, with appropriate principal axes ratio of about rj ~ 7.. 10, the sharp longitudinal 
resonance of nanoparticles can be tuned to the operational range of cladding modes 
resonances, thus making the TFBG sensor more sensitive to perturbations in the 
surrounding media. The sensor can even be coated with silver nanorods of various 
lengths, we still should expect that a fraction of these particle would have an ap¬ 
propriate aspect ratio and thus would contribute to an increase of the TFBG sensor 
sensitivity. We also should note that in the case of prolate particles with high aspect 
ratio the quasi-static approximation might not provide a satisfactory prediction, as 
was also the case with large spherical nanoparticles. 

7.3 The optimal parameters choice for coatings 
based on spherical nanoparticles 

Our goal is to design a nanoparticle coating with a sharp absorption resonance, 
located in the operational range of the TFBG sensor with location sensitive to the 
refractive index of surrounding media. From the previews section we know that silver 
or gold prolate nanoparticles might be an excellent choice, as the position of the 
plasmon resonance can be tuned at the exact location by choosing particles with the 
proper aspect ratio between its axes (77 ~ 7.. 10 ). 

In this section we discuss choice of parameters for coatings based on spherical 
nanoparticles. The sensor coated with spherical nanoparticles is schematically shown 
in Figure 7.14. 

Let us plot absorption efficiencies for various particle sizes and materials. As a first 
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Figure 7.14: A schematic representation of a TFBG sensor coated with spherical 
nanoparticles. 

step we will fix the the imaginary part of the particle permittivity, and compute Q a b.s 
as a function or real part of the relative complex refractive index '‘Rim] = = ^ n P arUcl ^ 
and the particle size parameter x = "y-R, as shown in Figures 7.15 and 7.16. 

Our goal here is to choose an optimal material and particle size. There are two 
main parameters to be considered: first we need a sharp prominent resonance, so 
that the real part of the effective refractive index can be significantly altered, and 
the second is the resonance steepness. We have a particular interest in the steepness 
of resonances. The steeper is the resonance the more its location is affected by the 
refractive index of the external medium, hence the higher is the sensor sensitivity. 

We can conclude that it is desirable to have particles made of material with a 
small imaginary part of the dielectric permittivity A[m] ~ 0.05 and choose a particle 
size parameter x such that 9ft [m] ~ 2 — 3.5. We note that the smaller is the 9ft [m] 
parameter, the steeper is the resonance, unfortunately for small 9ft[m] we can no 
longer obtain a prominent resonance. Once the real part of dielectric permittivity 
was chosen, we can study the effect of material absorption in more detail, as shown 
in Figure 7.17. We note that for 9f[m] < 0.05 the resonance amplitude starts to 
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lm(m) = 0.5 Q abs Q abs 



lm(m) = 0.3 Q abs Q abs 



Figure 7.15: The absorption efficiency as a function of the particle size parameter x 
and the real part of the relative complex refractive index 9ft [m] at various fixed values 
of 3f[m]. 


decrease, hence a material with small but noticeable loss is required. 
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Figure 7.16: The absorption efficiency as a function of the particle size parameter x 
and the real part of the relative complex refractive index 9?[m] at various fixed values 
of 3f[m]. 
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Figure 7.17: The absorption efficiency as a function of the particle size parameter x 
and the imaginary part of the relative complex refractive index 3f[m] at various fixed 
values of 9?[m]. 
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Considering the large value of the required relative refractive index m our choice is 
limited to only a few materials with high refractive indices, such us Zirconium Silicate 
and Titanium Dioxide. 

The refractive index of Titanium Dioxide is defined by equation (7.7) and is shown 
as a function of frequency in Figure 7.18. 


n = 5.913 + 


0.2441 
A 2 - 0.0803 


(7.7) 



Figure 7.18: Refractive index of titanium dioxide, TiO 2 . 

Considering the operational wavelength A = 1.55/im and the refractive index of 
the external medium ( n me di U m = 1-318 for water) we find that the relative refractive 
index of 7A0 2 particle immersed in water is m ~ 1.9. The horizontal line correspond¬ 
ing to m ~ 1.9, as shown in Figure 7.19, crosses several Mie resonances, hence a 
proper particle size can be chosen at one of such crossing points. Let us choose the 
particle size parameter x = 2.1 , corresponding to the second resonance. Then the 
particle radius can be found from the relation x = k 0 n m R and for the above mentioned 
parameter is R = 393 nm. 

We note that we only barely met the required conditions. With the available 
refractive index for particles we reached the very beginning of the resonance peak, 
as shown in Figure 7.19. A material with a higher refractive index in the IR band is 
highly desirable, unfortunately our choice is limited. 
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Figure 7.19: The optimal size of a particle made of Titanium Dioxide. 

The properties of the TFBG sensor coated with spherical nanoparticles can be 
modeled with equation (7.6) as previously. In the case of sparsely deposited spheri¬ 
cal particles the dipole-dipole radiative interaction between the particles can be ne¬ 
glected. Thus the effective dielectric permittivity function can be assumed to be an 
isotropic function e p = e<y Hence, once the particles absorption is known the re¬ 
fractive index can be computed with help of the Kramers-Kronig relation (5.1) and 
the Maxwell-Garnett model, used to compute the effective refractive index of the 
particles embedded in the surrounding medium. 

7.4 Conclusion 

In this chapter we presented results of our simulations of spherical and ellipsoidal 
nanoparticles made of various materials. We presented the method for increasing 
TFBG sensor sensitivity based on resonant coupling between the modes of TFBG 
sensor and the modes of particles deposited on the fibre surface. The particle are 
chosen to have a sharp resonance with its location sensitive to the external refractive 







Chapter 7. The optimal parameters for a nanoparticle-based coating 192 


index. 

It was shown that the metallic nanorod-like particles made of silver or gold look 
promising for the sensitivity-enhancing coating. The resonance of nanorod particles 
can be tuned exactly to the operational range of the sensor by choosing particles with 
77 ~ 7.. 10 aspect ratio between its principal axes. 

Alternatively, dielectric spherical particles with high refractive index can be cho¬ 
sen. We found that particles made of Titanium Dioxide with the diameter D ~ 800 nm 
look promising for the sensitivity enhancing coating. 

We should also note that the particle-based coatings should be sufficiently sparse, 
so that the interaction between particles, and hence the broadening of the resonance 
peak can be neglected. However, the sensitivity enhancement effect is proportional 
to the number of particles deposited on the fiber the surface. The effective refractive 
index of nanoparticle-based coatings can be computed with the Maxwell Garnett 
model (see Section 5.3), if the particles are assumed to be mutual independent and 
non-interacting. 

In the next chapter we present our experimental findings. 



Chapter 8 


Modification of the sensor surface 
with various types of nano-scale 
coatings 


In this chapter we present our experimental results on enhancing the TFBG sensor 
sensitivity by coating its surface with various types of nano-scale coatings. It is 
already known that the TFBG sensor performance can be improved by coating its 
surface with a nano-scale metal film [7, 15], thus coupling the TFBG resonances to 
a surface plasmon resonance (SPR) excited in the metal film. This technique was 
thoroughly investigated in [7, 15]. In the presented work we are mainly interested 
on a possible sensitivity improvement with the introduction of nanoparticle-based 
coatings. 

As was shown in Chapter 7, we propose using an external resonant system, such 
as small nanoparticles coupled to the TFBG resonances, to boost the sensor perfor¬ 
mance. It is well known that the position of a particle resonance is strongly affected 
by the surrounding medium, as shown in Figures 7.19 and 7.7, thus a thin coating 
layer created out of such particles, with optical properties strongly dependent on the 
refractive index of the external environment, can be engineered to boost the TFBG 
sensor sensitivity. For example, such effect was observed in [3] where the sensor was 
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coated with silver nanowires. A 3.5-fold increase in the TFBG sensor sensitivity was 
reported. In this context we have considered various types of particles. However, 
spherical and ellipsoidal particles are particularly interesting as the exact analytical 
solution exists for such cases. 

During the course of our research we investigated several coatings, such as metallic 
films of various morphologies and thicknesses, as shown in Figures 8.1 and coatings 
with nanoparticles of various shapes and materials. In particular we have tested coat¬ 
ings with nanocubes of various sizes (40 — 80 nm) synthesized from silver and gold 
metals, as shown in Figure 8.5, coatings with elongated gold and silver nanoparticles, 
as shown in Figures 8.2 and 8.3, and coatings with nanospheres, shown in Figure 8.4. 
The particles were deposited at various densities creating films of different morpholo¬ 
gies. The surface images were obtained by means of a scanning electronic microscope 
(SEM) and an atomic force microscope (AFM). 
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Figure 8.1: The SEM mages of the fibre surface coated with gold (top) and copper 
(bottom) films of different morphologies and thickness [1]. 




Figure 8.2: The AFM (a) and SEM (b) images of the fibre surface coated with silver 
nanowires [3]. 
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Figure 8.3: The SEM images of gold nanorod based coating 
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Figure 8.4: The SEM images of Tf0 2 spheres deposited on the fibre surface 



Figure 8.5: The AFM and SEM images of silver nanocubes (80nm) 
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8.1 Deposition and synthesis techniques 

In this section we briefly describe chemical methods of particles synthesis and 
deposition techniques that were used in the presented work. 

8.1.1 Synthesis of silver nanorods and nanowires 

Silver nanowires were chemically synthesized as described in [159]. The procedure 
results in highly crystalline nanowires with smooth surfaces. All the reagents used 
for synthesis were obtained from Sigma-Aldrich. All glassware was cleaned using 
aqua regia, rinsed in 18.2 MLl deionized water, and placed in an oven to dry prior to 
experimentation. 

A 50 mL round bottom flask containing 24.0 mL of anhydrous 99.8% ethylene 
glycol (EG) and a clean stir bar were placed in an oil bath set to 150° C and allowed 
to sit for 1 hour. Using a micropipette, 400 gL of 3 mM sodium sulfide dissolved 
in EG was added to the flask. Ten minutes later 6 mL of EG containing 0.12 g 
of dissolved poly(vinylpyrolidone) (PVP) with a molecular weight of 55000 amu was 
injected using a glass syringe immediately followed by the injection of 0.5 mL of 6 mM 
HC1. After an additional 5 minutes, 2.0 mL of 282 mM 99% + silver nitrate dissolved 
in EG was injected slowly using a glass syringe. Upon addition of the silver nitrate 
the solution immediately turned black and slowly became a transparent yellow, then 
changed to an ochre color as some plating in the flask occurred. The reaction was 
allowed to continue until the solution became white. 

The reaction was monitored by periodically taking small samples out of the re¬ 
action flask using a pasteur pipette and dispersing it in a quartz cuvette filled with 
95% ethanol for UV-visible spectroscopy. The reaction was quenched by placing the 
flask in an ice bath after the solution and had become fully white and turbid in ap¬ 
pearance. Purification of the nanowires was achieved by adding 20 mL of ethanol to 
the solution and centrifuging it at 13800 g for 20 minutes to remove the excess PVP, 
EG, and any reaction by-products. The supernatant was then discarded and the rods 
re-dispersed in ethanol by sonication. This process was then repeated several times 
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at 400 g to separate out the heavier wires from the solution. 

8.1.2 Synthesis of gold nanoparticles 

The gold nanoparticles were synthesized by the reaction of citrate reduction of 
chloroauric Acid (H[AuCl^\) in water, as was proposed by Turkevich [160, 161, 162], 
The reducing agent sodium citrate (Na 3 C 6 H 5 0 7 ) reduces gold ions to neutral gold 
atoms, when the solution becomes supersaturated gold nanoparticlcs start to form. 
Citrate ions act as both a reducing agent and a capping agent. The repulsion of 
the negatively charged citrate ions prevents the nanoparticles from forming aggre¬ 
gates [162]. As the result of this process spherical gold nanoparticles in the range 
from 10 to 20 nm in diameter suspended in water were created. 

8.1.3 Deposition of nanoparticles 

At the first step, the fibre surface was cleaned and prepared for chemical deposi¬ 
tion. The plastic jacket around the fibre was removed by soaking it into dichloromethane 
(CH 2 CI 2 ) and further cleaned to remove organic residue through the following mul¬ 
tistep approach. The uncoated area of the fibre was rinsed with methanol and sub¬ 
sequently treated with freshly prepared piranha solution (mixture of H 2 O : NH :i : 
H 0 O 2 = 5 : 1 : 1) at 70° C for 10 minutes to remove the remaining organic residues, 
and rinsed in 18.2 Mfl deionized water. 

The cleaned fibre was then submersed in a 1% (3-aminopropyl)trime-thoxysilane 
(APTMS, Aldrich, 97%, 281778) in methanol for 20 minutes in order to form a uniform 
a positively charged monolayer on the fibre surface. The APTMS molecules assembled 
on the glass by covalent bonding were exposed to hydroxyl sites (Si—OH) on the glass, 
thus forming a cross-linked self-assembled monolayer (SAM). The APTMS treatment 
replaces the hydroxyl groups (OH) adsorbed on the glass (fibre) substrate (S 2 O 2 ) with 
APTMS molecules forming a siloxane bond between the Si on one end of the APTMS 
molecules and an oxygen atom on the Si02 surface [163]. As a consequence, the amino 
group attached on the other end of the APTMS molecule is oriented away from the 
substrate. These amino groups on the APTMS molecules immobilise gold particles 
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onto the substrate because of the affinity of the amino group to the gold [163]. 

Finally, the APTMS modified fibre was rinsed with methanol and deionized wa¬ 
ter followed by drying with a flow of gas. After drying, the modified fibre was 
submerged into a freshly prepared colloidal solution of gold nanoparticles and left 
for 24 hours. The results of nanoparticles deposition are shown in Figures 8.2, 8.3 
and 8.5,. 

8.1.4 Electroless metal coating 

The fact that a cylindrically symmetric fibre is used makes the deposition of an 
uniform nanometer-scale him quite challenging. Previously in our group the time 
calibrated two-step sputtering process was used in order to deposit gold films [7, 164], 
but we faced some reproducibility issues that could only come from coating errors. 
In an effort to improve this situation we had investigated electroless plating [165] 
because of its potential for highly uniform metal coatings on non-planar surfaces, its 
low cost and simplicity of implementation (with a potential for batch production), and 
its relatively low deposition rate, which facilitates the control of the process duration. 
The electroless gold and silver coatings can be deposited reliably and accurately for 
the fabrication of near infrared TFBG-SPR sensors, but most importantly the TFBG 
itself can be used to monitor the deposition process and to stop it at the optimum 
him thickness, for example for SPR excitation of metal coated hbres in water. While 
coating conventional FBGs with copper and nickel had been described earlier for the 
purpose of making the FBG response athermal, the exact thickness and the uniformity 
of the films was not critical for those applications [165]. Further studies of thermal 
stress between such metal coatings and optical hbres were performed in [166], where it 
was shown that the metal coating of FBG sensors can be used not only as a protective 
layer but also as a temperature sensitivity booster. 

The method we were using to coat the TFBG sensor relies on the attachment of 
gold nanoparticles on a suitably prepared bare hbre surface, followed by the electroless 
plating process, based on the reduction of metallic ions from a solution to a solid 
surface without applying an electrical potential [167, 168]. 
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To initiate the reaction the fibre surface was precoated with gold nanoparticles, 
with the method described in the previous section. The electroless plating is based 
on the process of negative charge production on the surface of gold nanoparticles 
through the oxidation reaction assisted by a reducing agent. The gold nanoparticles 
immobilized on the glass surface constitute excellent sites for the reduction of gold or 
silver ions by the hydroxyl amine solution. Thus, the nanoparticles attached to the 
fibre surface play a catalytic role during the process of electroless plating. 

The modified with gold nanoparticles fibre was next dipped into the plating bath. 
For the silver coating the plating solution consisted of silver nitrate ( AgNO 3 , 0.01 M), 
ammonium nitrate ( NH 4 NO 3 , 8.96M), acetic acid (CH 3 COOH , 2.24M), and hy¬ 
drazine hydrate (. H 2 NH 2 .H 20 , 0.4 M) at a volumetric ratio of 1 : 1 : 1 : 2, re¬ 
spectively. For the gold plating the solution consisted of 0.01% chloroauric acid 
(HA 11 CI 4 . 3 H 2 O) and the reducing agent hydroxylamine hydrochloride ( NH 2 OHHCl , 
0.4 M ) in a 1 : 1 volumetric ratio. 

The electroless plating occurs on the surfaces of already immobilized gold nanopar¬ 
ticles, eventually merging into a continuous thin metallic him, as shown in Figure 8.1. 
The him thickness can be controlled in real-time with high precision [1], The him 
morphology can also be controlled by varying the deposition time, the solution com¬ 
position, or the type and size of nanoparticles. 

As a separate check we further calibrated the monitoring process by performing 
Atomic Force Measurements (AFM) on sample films to correlate the TFBG response 
with the physical thickness of the films. The AFM measurements provided evidence 
that the electroless-deposited gold layer retained some significant granularity. During 
the deposition, we monitored the wavelengths and amplitudes of several resonances 
in the transmission spectrum. 
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8.2 The optimal metal coating for Surface 
Plasmon Resonance (SPR) excitation 

The sensitivity of TFBG sensor can be enhanced by coating its surface with a 
thin metal film, of a proper thickness, thus allowing the energy coupling between the 
backward-propagating cladding modes and the collective electron oscillation on the 
metal surface, the so-called surface plasmon resonances (SPR) [13]. The sensitivity 
enhancement of the TFBG-SPR sensor was previously reported in [7, 15]. Here we 
focus on optimization of the coating process. 

In order to optimize the SPR sensitivity and resolution, the metal coating thick¬ 
ness and uniformity must be controlled very accurately around the fibre circumference, 
for thicknesses of the order of a few tens of nm. For coatings that are too thin, the 
SPR resonances broaden and weaken. On the other hand, when the thickness is too 
large, light cannot tunnel across the metal to excite the plasmon polariton on the 
outer surface. 

Here we present a novel method for the optimal metal coating of TFBG sensors 
by electroless plating, thus creating a Surface Plasmon Resonance (SPR) finely tuned 
into the operational range of the TFBG sensor [1], 

The electroless plating of gold or silver metals from solution is a simple and efficient 
way to deposit a conformal coating of the required thickness and uniformity for fibre 
SPR applications, and furthermore, the deposition process can be controlled in real¬ 
time by means of a TFBG sensor and stopped at the optimum film thickness for SPR 
excitation. Therefore, a precise control of the plating process (temperature, chemical 
concentrations) is not necessary and reproducible coatings with thicknesses on the 
order of 50 nm can be obtained regardless of the speed of the plating process. 

A quantitative measure of the quality of the SPR response of the TFBG in the 
plating solution is provided by tracking the polarization-dependent loss (PDL) in real 
time and used as a criterion to end the plating. When the metal coating has a thick¬ 
ness for which a plasmon wave can be excited by cladding modes, the PDL spectrum 
acquires a deep notch that reveals the SPR location [15]. Using the polarization- 
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based optical sensing method, presented in Chapter 4, we can identify the precise 
moment of the plating process at which the SPR signature is observed. When the 
spectral notch is maximized, the thickness is optimal for SPR operation. This is best 
observed in the density plot of Figure 8.6, where the envelopes of the PDL spectra 
are represented as a colour scale in the horizontal direction at different moments of 
the plating process, the vertical direction. 



Figure 8.6: (a) The envelope of PDL spectra, taken continuously along the course of 
gold him deposition, and cross sections centered at the point of the deepest notch: 
wavelength = 1542 nm (b) and time = 7 min (c) 


Indeed, it can be clearly seen that after only approximately 7 minutes of gold 
him deposition a deep notch in the PDL envelope occurs (darkest blue colour in 
Figure 8.6(a)). The individual PDL spectrum corresponding to this plating time is 
shown as Figure 8.6(b), while Figure 8.6(c) extracts the time evolution of the PDL 
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value at the wavelength where the SPR is observed. The optimum point is clearly 
seen in the latter figure, hence the deposition process can be interrupted as soon as 
a local PDL envelope minimum is detected, regardless of the plating rate. 

We should note that the SPR signature appears suddenly, hence the monitoring 
process is necessary. Figure 8.6(c) further shows that the PDL envelope evolution 
slows down gradually for longer plating times. These effects are due to both the 
self-termination of the plating process and to the fact that as the metal layer be¬ 
comes thicker it eventually shields the light from the cladding modes from the outside 
medium so that further thickness growth is not detected. 

We conclude that using TFBG’s, we have presented a new method to monitor the 
growth of the plated film in situ and in real time. Most importantly however, we 
showed that by monitoring a simple parameter in the polarization-dependent Loss of 
the TFBG during plating, the optimum film thickness for SPR operation can be found 
with great accuracy, regardless of the plating rate, ft is also likely that other metals 
could be plated and monitored in similar fashion. The method is inherently limited 
in thickness since the evanescent field of the cladding modes must tunnel across to 
the outer surface of the metal film in order to detect changes in thickness. When the 
thickness exceeds the penetration depth of the light, the thickness growth appears to 
saturate. Obviously for the application of fibre SPR, the thickness needed must be 
smaller than the penetration depth, so this does not present an actual limitation in 
this case. 


8.3 Sensitivity enhancement with a nanoparticle 
based coating 

As we discussed in the previous section the sensitivity of a TFBG sensor can 
be increased by coating its surface with a metal film of the appropriate thickness, 
thus creating conditions for SPR excitation. However, such improvement occur only 
over a narrow spectral range for which the phase matching condition between fibres 
cladding modes and the SPR wavevector is satisfied. In this section we present a new 
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technique for increasing the sensor sensitivity. Instead of a resonance in continuous 
metallic him we proposing to use resonances of nanoparticles deposited on the sensor 
surface [3]. 

We showed in [3] that by sparsely coating the TFBG sensor with randomly oriented 
silver nanowires, a large number of fibre modes with different wavevectors and electric 
held distribution can resonate with nanowires, exciting localized surface plasmon 
resonances (LSPR). LSPR are collective oscillations of conducting electrons in the 
nanoparticles giving rise to strongly enhanced and highly localized electromagnetic 
holds, which can be utilized in various sensing devices [7]. 

8.3.1 Coating with silver nanowires 



SEM MAG: 10.00 kx SEM HV: 20.00 kV 1. -1 -1 — 1 

Del BSE WD: 7.251 mm \ 

View field: 15.00 pm Date<m/d/y): 01/10/12 1 ^ 


Figure 8.7: The SEM (b) images of the hbre surface coated with silver nanowires [3]. 

The sensor was coated with chemically synthesized silver nanowires ~ 100 nm 
in diameter and several micrometres in length, as shown in Figure 8.2. A UV-vis- 
NIR absorption spectrum was collected on a similar sample deposited on a hat glass 
slide. The spectrum shown in Figure 8.8 contains two main features: a peak at ~ 
380 — 400 nm corresponding to the excitation of transverse plasmon resonances along 
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the short axis and a broad peak in the NIR region corresponding to the excitation of 
longitudinal plasmon resonances along the long axis. 



Wavelength (nm) 

Figure 8.8: UV-vis-NIR absorption spectrum of synthesized nanowire coating (de¬ 
posited on a flat glass substrate). The insets illustrate the relative polarization of the 
Plasmon oscillations that give rise to the absorption. 

The spectral region of interest for TFBG (A E [1525,1590] nm) is highlighted in 
grey on the figure. It is seen that a strong extinction signal for silver nanowires exists 
at these frequencies. However, since the nanowires were deposited on the glass surfaces 
(fibre or flat substrate) using the self-assembly approach [169] (made possible by their 
partially negative charge [170]), their orientation was not controlled (Figure 8.2). 
Therefore the excitation of the longitudinal plasmons in the 1550 — 1600 nm band by 
polarized light sources can only occur for those wires that are parallel to the incident 
light polarization (i.e. about half on average). It appears impossible to excite the 
short (radial) plasmons of these nanowires at the same wavelengths. 

At the first step of the characterization process, the TFBG sensor was tested to 
determine its refractometric operational range before and after the nanowire coating 
deposition. Mixtures of ethylene glycol (EG) and water were used with different 
volume ratios, varying from pure water to pure EG, allowing to cover a range of 
An = 3.81 x 10~ 2 for the refractive index of the surrounding medium. A custom- 
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made Teflon cell was used to hold the fibre still and straight while being submerged in 
varying solutions of ethylene glycol (EG) and water while transmission measurements 
were carried out with an optical vector analyser (OVA). Figure 8.9 shows the average 
insertion loss of the coated and uncoated sensor for various intermediate values of the 
surrounding refractive index. 
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Figure 8.9: The TFBG spectrum evolution before (a) and after (b) deposition of 
nanowires, for several values of the refractive index of the solution. The concentration 
of the Ethylene Glycol in water goes from 0%, 25%, 50%, 75%, 100% from the top 
to the bottom of the figure. The corresponding total refractive index change is An = 
3.81 x 10 -2 . 


In those spectra, each downward peak corresponds to loss of light in the core, due 
to coupling to cladding modes. A discontinuity in the amplitudes of the resonances 
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is observed on all spectra and its spectral position shifts towards longer wavelengths 
as the surrounding index increases. This discontinuity corresponds to the cut-off con¬ 
dition where the cladding modes become leaky because their effective index becomes 
larger than the surrounding index [11]. The behavior of both groups of spectra is 
similar, apart from the fact that the amplitudes of the resonances are always smaller 
for the coated grating. This is consistent with the assumption that plasmons can be 
excited in the nanowires in the 1525 — 1590 nm range, thereby inducing additional 
loss in the cladding modes, which results in a decrease of the coupling strength (and 
of the resonance amplitude). The observation of the shift in the cut-off wavelength of 
the average insertion loss provides a relatively rough (but absolute) estimate of the 
refractive index of the solution. 

For finer measurements (of small index change increments), observations on a sin¬ 
gle resonance provides better accuracy. The refractometric sensitivity measurements 
were performed on the TFBGs before and after the deposition of silver nanowires 
by monitoring the response following small changes of the relative EG concentration. 
The initial measurements were taken using a well mixed 1 : 1 mixture of EG and H20 
prepared in a large 200 mL beaker and transferred to the 5 ml Teflon cell. Subsequent 
measurements were taken with solutions obtained by adding 10 /il of water to the 
cell and continuous stirring. The mechanical stirring was stopped prior to each mea¬ 
surement to eliminate potentially detrimental fluid motions and acoustic disturbances 
around the fibre that could result in false readings. After each dilution, the refractive 
index of the solution was determined from Lorentz-Lorenz mixing equation 5.24 [171]: 


n — 1 717 — 1 

- =P'-TT T+P2 2 


n 2 2 -l 


( 8 . 1 ) 


n 2 + 2 J 1 ri{ + 2 r ‘ri2 + 2 
Here, n represents the refractive index of the mixture, ri\ and ri 2 are refractive indices 
of pure components (water and EG, 1.318 and 1.394, respectively at wavelengths 
near 1550 nm [172]), and p\ = Vl/V, P 2 = V2/V are the volume fractions of the 
components. 

Figure 8.10 (a and b) depicts the response of an individual resonance to small index 
variations, before and after the sensor surface was modified with silver nanowires, and 
for the two singular values of the transfer function. 
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Figure 8.10: Singular values (a,b) and polarization-dependent loss parameter (c,d) 
(linear scale) changes due to a small refractive index change of An = 3.77 x 10“ 4 , 
before (a,c) and after (b,d) deposition, corresponding to a single resonance taken at 
A = 1555.7 nm . 


We note that the most obvious effect resulting from the nanowire coating is the 
increase in the splitting between the two singular values, from ~ 30 pm to ~ 100 pm 
(comparing the central wavelengths of Ai and A 2 for n\ or 712 ). This increase should be 
associated with an increase in the polarization dependence of the boundary condition 
for the cladding modes, as expected for metal particles. A similar effect had been 
observed previously by our group for gratings covered by a sparse layer of spherical 
copper nanoparticles deposited by a pulsed chemical vapor depositio (CVD) tech¬ 
nique [173]. 

It was also observed that the singular value spectra shifted to shorter wavelengths 
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with decreasing refractive index. However, this shift is somewhat difficult to quantify 
precisely because the resonances are rather broad and flat bottomed. Instead of 
isolating one of the singular value to quantify the shifts (and thereby throwing away 
half the data, corresponding to the other singular value) we will follow zeros in PDL 
spectra, corresponding to resonances of interest, as we discussed in Chapter 4. We are 
basically measure differential spectrum between polarization states, which provides 
extremely sharp peaks (Figure 8.10c and 8.10d) that can be tracked down with much 
greater accuracy than the insertion loss or singular value resonance position. 

Conducting subsequent measurements for fine variations in the refractive index 
obtained by the step-wise addition of 10 /if quantities of water to the solution we 
determine the refractometric sensitivity for several peaks. Each addition of water 
caused a small decrement of An ~ 7.5 x 10 -5 in the refractive index of the surround¬ 
ing medium. Those measurements can be carried out on any of the resonances. The 
results of these detailed measurements are shown in Figure 8.11, where a set of 5 
different individual resonances were measured in order to determine if the spectral 
position of a resonance had an impact on its sensitivity, shown in Figure 8.11b. 

Comparing the slopes in Figure 8.11b for the resonances of the coated and un¬ 
coated TFBG, we obtain an increase in sensitivity from 53 run RIU -1 for the bare 
uncoated TFBG sensor to 185 rim RIU -1 for the sensor coated with silver nanowires 
(RIU stands for “Refractive index unit”). Although the reported sensitivity for metal 
coated SPR TFBG sensors is reported to have a higher value, of about 555 run 
RIU~ l [174], the operational range of SPR sensors is limited to a narrow range of 
few nanometers where the resonant phenomena of SPR excitation is observed. The 
proposed sensor with a coating of nanowires has an advantage over a SPR sensor as 
the increase in sensitivity is observed over the whole operational range of the sensor, 
as shown in Figure 8.11a. 
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Figure 8.11: (a) Wavelength shift AA of several individual resonances of the TFBG 
sensor due to the surrounding refractive index change An. (b) The sensetivity before 
(open marks) and after (closed marks) deposition [3]. 
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8.3.2 Discussion 

The observed enhancement in sensitivity was investigated by taking a close look 
at the electromagnetic field interaction with the silver nanowires. The device is ca¬ 
pable of sensing the external medium by means of the evanescent field of cladding 
modes, leaking outside the fibre cladding, as shown in Figure 3.46. We note that the 
evanescent field decays slowly enough to extend into the thin layer of silver nanowires, 
thus conditions for coupling between the evanescent field incident from the fibre and 
the nanoparticles on its surface are created. Perturbations of the evanescent field 
modify how the grating couples the light from the core to each cladding mode and 
hence the resonances observed in the transmission spectrum. Knowledge of the fi¬ 
bre geometry and refractive indices, as well as of the grating period, tilt angle and 
modulation amplitude, allow us to uniquely assign specific cladding modes to the 
resonances observed experimentally. The electric field distribution corresponding to 
two closely positioned TM-like and TE-like modes of the same azimuthal family at 
A = 1548.1 nm and A = 1552.1 nrri are shown in Figure 8.12. 

a) b) 



Figure 8.12: Vectorial E field structure for two modes with almost identical prop¬ 
agation constants (hence resonance wavelengths), but different polarization states 
((a)TM-like mode, (b)TE-like mode). Silver nanowires are shown schematically on 
top. 
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As discussed in Section 7.2.3, the coatings with cylindrical particles give rise to 
very prominent and sharp longitudinal resonances, which shift toward the infrared 
band with an increase of particle elongation. The transversely polarized electric 
held excites a series of blue-shifted transverse resonances in the case of elongated 
nanoparticles. The particles should be chosen in such a way that the resonances are 
positioned far away from the operational range of the sensor, and hence the dielectric 
permittivity e p is a relatively constant function in the operational range. However, 
the longitudinal resonances of the particles, excited by the tangentially polarized 
electric held E are positioned exactly in the range where the sensor operates, so 
that the energy coupling can occur between the TE-likc modes of the sensor and local 
resonances of the particles, as shown in Figure 8.12. The imaginary part of the 
dielectric permittivity function has a resonance in the operational range. The position 
of the resonance is defined, among other parameter, by the external refractive index 
of the surrounding medium. Therefore the differential sensitivity with the TE-like 
and TM-like modes is possible, due to the fact that the modes are affected differently 
by changes in the external refractive index. 

The TFBG couples core-guided light to two different families of cladding modes 
according to the polarization of the input light relative to the orientation of the tilt 
plane, as shown in Figure 3.47. To be precise, when the input core-guided light 
is polarized linearly in the plane of the tilt (corresponding to P-polarized light), 
the cladding modes that are excited by the TFBG belong to the TM-like family of 
modes and have their electric held polarized predominantly in the radial direction at 
the cladding boundary. On the other hand, the S-polarized light (again relative to 
the plane of tilt), and the grating couples the core light to TE-like cladding modes 
whose electric held at the cladding boundary is mostly azimuthal (i.e. tangential). 
The calculated wavelength spacing and refractometric properties of the two modes 
in each pair reveal that they can be associated with the observed singular values of 
the insertion loss resonances measured above. Therefore, we can explain the observed 
splitting of the singular value pair by noting how the electric held of the correspond¬ 
ing mode probes the nanowires: since the TM-like modes have electric helds that 
are predominantly radial at the cladding boundary while TE-likc modes are mostly 
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tangential, TM-like modes cannot excite the “long axis” plasmons of the nanowires 
while the TE-like modes can. 

We also note that these results have been obtained in a relatively stable tem¬ 
perature environment, even though no active temperature control was used. In 
was reported [11] that the underlying TFBG platform can be made temperature- 
independent over several tens of degrees by referencing all wavelengths to the Bragg 
resonance, bounded the fibre’s core and completely insensitive to changes outside the 
fibre. With this referencing scheme the only impact of temperature on the positions of 
the resonances arises from the change in the refractive index of the nanoparticles and 
of the medium that surrounds them. The temperature effect was studied with appli¬ 
cation to TFBGs coated with uniform gold films and it was found to be insignificant 
(about 10 pm /°C [175]) compared to the shifts observed here. 

8.3.3 Conclusion 

The sensitivity of a TFBG refractometer can be increased at least 3.5-fold by the 
addition of a sparse coating of silver nanowires. The coating is created by simple liquid 
phase self-assembly process that does not require special deposition tools, unlike other 
methods used to fabricate plasmon assisted devices (e-beam lithography, chemical 
vapour deposition [7, 176, 86] or real time monitoring of the deposition process [1] to 
achieve very specific layer thicknesses). 

The high sensitivity of the sensor is explained by the fact that the resonances of 
a l(cm) long grating have Q values in excess of 15000 at near infrared wavelengths, 
due to full widths of 100(pm). As a result, the positions of these resonances can be 
followed accurately, even for shifts of the order of several picometres (reproducibility 
better than 3 pm was recently demonstrated experimentally with a similar measuring 
system [15]). 

Associated with our 185 nm RIU -1 sensitivity, a 3(pm) resonance position accu¬ 
racy yields a minimum detection level of 1.6 x 10 5 RIU. 

The second advantage is based on the fact that different resonances have simi¬ 
lar sensitivity and any of them can be used for refractive index sensing, as opposed 
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to other kinds of SPR sensor, where only one or a few resonances satisfying the 
phase matching condition have high sensitivity [15]. The high Q value of these reso¬ 
nances is a key factor in taking advantage of plasmonic effects because of the inherent 
trade-off between strong plasmon excitation (and enhanced sensitivity) and increased 
loss, which tends to broaden resonances: the TFBG configuration allows the device 
designer a great amount of control to adjust the amount of overlap between the 
evanescent held and the metal particles and also in the relative orientation of the 
light polarization and particle geometry. 



Chapter 9 


Conclusion 


The primary focus of this book has been set on theoretical analysis of the TFBG 
sensor and proposing methods of enhancing its sensitivity by surface modification. 

The TFBG sensor is not a trivial object. The sensor has a complex polarization- 
dependent response dependent on more than a thousand interacting modes. To ad¬ 
dress the challenge of the theoretical analysis, we developed a highly efficient and fast 
numerical solver, capable of computing a TFBG spectra in approximately 3 minutes 
for a given state of incident light polarization. 

Theoretical analysis 

The solver we developed consists of two major parts: the mode solver and the 
numerical integrator, which solve a system of coupled differential equations. 

The mode solver is the key element, as it is required to compute more than a 
thousand modes at various wavelengths. We developed a simple yet efficient and fast 
full-vectorial complex mode solver, capable of handling waveguides of an arbitrary 
complex refractive index profile. First, the system of Maxwell’s equations was reduced 
to only two coupled ordinary differential equations for the electric field. Next, the 
equations were transformed into a system of algebraic equations with the help of a 
finite difference method. Finally the problem was reduced to the problem of finding 
eigenvalues and eigenvectors of a five-diagonal sparse matrix. The eigenvalue problem 
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was effectively solved with the standard iterative method. 

At the next step the modes obtained for the non-perturbed case were used as 
the basis functions to solve the problem of tilted grating inscribed along the optical 
axis. The problem of small perturbations along the optical axis was handled with 
the coupled mode theory. As a result, the problem was reduced to a set of coupled 
differential equations, representing the energy transfer between the modes. Due to 
the grating tilt, the energy transfer depends on the core mode polarization. The 
polarization-dependent effects were thoroughly analysed in this work. Interesting 
insights into the properties of the electric field at the TFBG sensor boundary were 
gained. 

The results of the simulation were shown to be in a good accordance with the 
experimental measurements. 

Polarization-based experimental measurements 

Along with the theoretical study, we conducted an extensive experimental study of 
the optical properties of the TFBG sensor. We based our study on the Jones matrix 
and Stokes vectors data made available by means of Optical Vector Analyser (OVA). 

We proposed two alternative approaches towards the TFBG sensor data analysis. 
The first method was based on tracking the grating transmission of two orthogonal 
states of linear polarized light that were extracted from the measured Jones matrix or 
Stokes vectors of the TFBG transmission spectra. The second method was based on 
the measurements taken along the system principle axes and polarization-dependent 
loss (PDL) parameter, also calculated from measured data. It was shown that the 
frequent crossing of the Jones matrix eigenvalues as a function of wavelength leads 
to a non-physical interchange of the calculated principal axes. A method to remove 
this unwanted mathematical artefact and to restore the order of the system eigen¬ 
values and the corresponding principal axes was provided. A comparison of the two 
approaches revealed that the PDL method provides a smaller standard deviation and 
therefore a lower limit of detection in refractometric sensing. 

Furthermore, the polarization analysis of the measured spectra allows for the iden- 
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tification of the principal states of polarization of the sensor system and consequently 
for the calculation of the transmission spectrum for any incident polarization state. 
The stability of the orientation of the system principal axes was also investigated as 
a function of wavelength. A small oscillation in the orientation of the principal axes 
as a function of wavelength was observed. 

Sensitivity enhancement 

In the presented work we investigate the sensor response to various types of nano- 
scale him coatings. We have proposed the method of resonant coupling between the 
high quality factor resonances of the TFBG structure and the local resonances of 
nanoparticles deposited on the sensor surface. The problem was investigated theo¬ 
retically for spherical and elliptical particles made of various materials. The optimal 
parameters for particles were determined. A 3.5-fold increase in the TFBG sen¬ 
sor sensitivity was observed experimentally by coating the sensor surface with silver 
nanowires. 


This doctoral project resulted in the following publications: 

1 . A. Bialiayeu, C. Caucheteur, N. Ahamad, A. Ianoul, and J. Albert, ”Self-optimized 
metal coatings for fiber plasmonics by electroless deposition,” Optics express 19, 
18742-18753 (2011). 

In this paper we presented a novel method to prepare optimized metal coatings for 
infrared Surface Plasmon Resonance sensors by electroless plating. We show that 
Tilted Fiber Bragg grating sensors can be used to monitor in real-time the growth 
of gold nanohlms up to 70 nm in thickness and to stop the deposition of the gold 
at a thickness that maximizes the SPR (near 55 nm for sensors operating in the 
near infrared at wavelengths around 1550 nm). The deposited films are highly 
uniform around the fiber circumference and in spite of some nanoscale roughness 
(RMS surface roughness of 5.17 nm) the underlying gratings show high quality 
SPR responses in water. 
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2. A. Bialiayeu, A. Bottomley, D. Prezgot, A. Ianoul, and J. Albert, “Plasmon- 
enhanced refractometry using silver nanowire coatings on tilted fibre Bragg grat¬ 
ings,” Nanotechnology 23, 444012 (2012). 

In this paper we presented a novel technique for increasing the sensitivity of tilted 
fibre Bragg grating (TFBG) based refractometers. The TFBG sensor was coated 
with chemically synthesized silver nanowires ~ 100 nm in diameter and several mi¬ 
crometres in length. A 3.5-fold increase in sensor sensitivity was obtained relative 
to the uncoated TFBG sensor. This increase was associated with the excitation of 
surface plasmons by orthogonally polarized fibre cladding modes at wavelengths 
near 1.5 pm. Refractometric information was extracted from the sensor via the 
strong polarization-dependence of the grating resonances using a Jones matrix 
analysis of the transmission spectrum of the fibre. 

3. W. Zhou, D. J. Mandia, M. B. Griffiths, A. Bialiayeu, Y. Zhang, P. G. Gordon, S. 
T. Barry, and J. Albert, “Polarization-dependent properties of the cladding modes 
of a single mode fiber covered with gold nanoparticles,” Optics express 21, 245-255 
(2013). 

In this paper the properties of the high order cladding modes of a standard optical 
fiber were measured in real-time during the deposition of gold nanoparticle layers 
by chemical vapor deposition. A correlation between the transmission spectra of 
the 10° TFBG and the optical properties of gold particle coatings was established. 

My contribution to this work allowed to explain the observed effects by the nu¬ 
merical FDTD modelling of the gold particle coating layer, as described in Sec¬ 
tion 5.3.2. 

4. A. Bialiayeu, A. Ianoul, and J. Albert, “Engineering a resonant nanocoating for 
an optical refractive index sensor,” in “Electronic, photonic, plasmonic, phononic 
and magnetic properties of nanomaterials,” , vol. 1590 (AIP Publishing, 2014), 
vol. 1590, pp. 68-70. 

In this work we proposed an idea to boost the performance of refractive index 
sensor based on the tilted fiber Bragg grating structure by resonant coupling of 
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small spherical nanoparticles to the TFBG resonances. The optimal choice of 
nanoparticle parameters was discussed. 

5. A. Bialiayeu, J. Albert, A. Ianoul, A. Bottomlcy, and D. Prezgot, “Silver nanowire 
coated tilted fibre Bragg gratings,” in “Bragg Gratings, Photo-sensitivity, and 
Poling in Glass Waveguides,” (Optical Society of America, 2012), pp. BW2E1. 

In this work we reported a 3.5-fold increase in sensitivity of TFBG based refrac- 
tometers by coating the sensor surface with silver nanowires. A strong polarization 
dependence of the grating resonances was observed and analyzed. 
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MatLab Code. The full vectorial complex mode solver. 
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The following code had been developed to obtain the exact full vectorial solution 
to the problem of cylindrical and slab waveguides of an arbitrary refractive index 
profile, including lossy waveguides with a non-zero imaginary part of the refractive 
index. 

The main routine. 


1 % r in [—R, R] 

2 

3 clear all 

4 clc 

5 elf 

6 

7 %============== Input ===================================== 

8 lam =1.6; % wavelength um 

9 ko = 2*pi/lam; % free space wave vector 

10 

11 % layer vectors 

12 %RL = [4.1, 62.5, 70]; 

13 %nL = [1.4492, 1.444, 1.3556]; 

14 %nL = [1.4492+ 0.001, 1.444 + 0.001, 1.33]; 

15 

16 R = 1.5; 

17 % nl = 1.4492; n2 = 1.3556; % refractive index 

18 nl = 3; n2 = 1.3556; % refractive index 

19 RL = [R, R + 3]; 

20 nL = [nl, n2] ; 

21 

22 m=l; 

23 Neig =40; % number of eigenvalues to compute 

24 N = 10000; 

25 

26 %============== Potential ================================= 

27 

28 h = RL(end)/(N—0.5); 

29 r = linspace(0.5*h, RL(end), N) ; 

30 r = [— fliplr(r), r] . ' ; 

31 r = r — 0.1 *h; 

32 

33 ni = f.construct Jii2 (RL, nL, N) ; 

34 ni = [fliplr(ni), ni].'; 

35 

36 e — ni . ~ 2; 

37 %============== Solve =============================== 

38 

39 %[Ls, U] = f_FD_Construct_Matrix_Cyl_Exact_Vect (r, e, ko, m) ; 

40 [Ls, U] = f_FD_Construct_Matrix_Cyl_WGA_Vect (r, e, ko, m) ; 

41 %[Ls, U] = f_FD_Construct_Matrix_Cyl_Exact_Vect (r, e, ko, m) ; 

42 [V, Eig] = f_FD_eig(Ls, Neig, U) ; 
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43 

44 %============== Results =================================== 

45 Neff = sqrt(Eig./ko~2); 

46 U = sqrt (U/ko~2) ; 

47 

48 f_Plot_vec (m, r, U, Neff, V, 3, 1, 3.2) 

49 % f _P lot_Scal_Single (12, 3, —7, 7, 2.5, 3.1, m, r, sqrt(U)./ko, sqrt (Eig) ./ko, V, h) 


The matrix assembly for the weakly guided approximation of the scalar 
problem. 


1 function [Ls, u] = f_FD_Construct_Matrix_Cyl_WGA_Scalar (r, e, ko, m) 

2 % My working code 

3 % r — grid vector 

4 % ki = k_o n 

5 % m — if zero we get symmetric function: d_r u(r) = 0 

6 

7 

8 h = r(3) — r(2); % the step size 

9 N = length(r); % number of elements r 

10 

11 

12 %================== Construct d_x matrix =========================== 

13 B = zeros (N, 3) ; 

14 B (:, 1) = -1; 

15 B (:, 3) = 1; 

16 D = 1/(2*h)*spdiags(B, [ —1,0,1],N, N) ; 

17 

18 %================== Construct d_x ~2 matrix ========================= 

19 B = ones (N, 3) ; 

20 B (:, 2) = —2; 

21 D2 = 1/(h"2)*spdiags(B, [ — 1,0,1],N, N); 

22 

23 

24 %========================= Construct matrix ======================== 

25 % (d_r~2 + 1/r d_r + [k~2 - nT2/r~2 ] ) u (r) = b~2 u(r) 

26 

27 u = e*ko~2 — m~2 . /r. ~2; 

28 U = spdiags(u,0,N,N); 

29 

30 % sparse matrices 

31 Ri = spdiags(1./r,0,N,N); 

32 

33 Ls = D2 + Ri*D + U; 
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The matrix assembly for the weakly guided approximation of the 
vectorial problem. 


1 function [Ls, u] = f_FD_Construct_Matrix_Cyl_WGA_Vect (r, e, ko, m) 

2 

3 h = r(3) — r(2); % the step size 

4 N = length(r); % number of elements r 

5 

6 %================== Construct d_x matrix ========================= 

7 B = zeros (N,3); 

8 B(:,l)= -1; 

9 B (:, 3) = 1; 


10 D = 1/ (2*h) *spdiags (B, [ —1,0,1],N,N); 

11 

12 %================== Construct d_x ~ 2 matrix ==== 

13 B = ones (N, 3) ; 

14 B (:, 2) = —2; 

15 D2 = 1/(h~2) *spdiags (B, [ —1,0,1], N,N); 

16 

17 %================== Construct diagonal matrices 

18 % vectors 

19 u = e*ko~2 — rrT2 . /r. ~2; 

20 

21 ul = —(l+m~2) ./r."2 + e*ko~2; 

22 u2 = (m*2)./r."2; 

23 

24 % sparse matrices 

25 Ri = spdiags(1./r,0,N,N); 

26 Ul = spdiags(ul,0,N,N); 

27 U2 = spdiags (u2, 0,N,N) ; 

28 

29 L = D2 + Ri*D; 

30 Ls = [L+Ul, U2; U2, L+Ul]; 


The matrix assembly for the exact solution of waveguides with 
cylindrical symmetry. 


1 function [Ls, u] = f _FD_Const ruct_Matrix_Cyl_Exact_Vect (r, e, ko, m) 

2 % puled r under E_r, Need it as the dispersion curves are more realistic !!! 

3 

4 h = r(3) — r(2); % the step size 

5 N = length(r); % number of elements r 

6 

7 %============== 

8 B = zeros(N,3); 

9 B (:, 1) = -1; 


Construct d_x matrix 
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10 

B (: 

, 3)= 1; 

11 

D = 

1/(2*h)*spdiags(B, [ —1,0,1],N,N); 

12 



13 

%— 

- Construct d_x 2 i 

14 

B = 

ones (N, 3) ; 

15 

B (: 

, 2) = -2; 

16 

D2 

= 1/(h~2)*spdiags(B,[ —1,0,1],N,N) 

17 



18 

%== 

==================== The operator 

19 

ri 

= 1. /r; 

20 

ri2 

= ri."2; 

21 



22 

u = 

e*ko' N 2 — nT2*ri2; 

23 

U = 

spdiags (u, 0,N,N); 

24 

Ri 

= spdiags(ri,0,N,N); 

25 



26 

%— 

- TE—like modes 

27 

L2 

= D2 - Ri*D + U; 

28 



29 

%- 

- TM—like modes 

30 

ei 

= 1./e; 

31 

de 

= D*e; 

32 

dine = de.*ei; 

33 

P = 

ri + dine; 

34 

P = 

spdiags(p,0,N,N) ; 

35 

LI 

= D2 - P*D + U; 

36 



37 

% — 

- off diagonal terms 

38 

a = 

2 *m*ri2; 

39 

ul 

= e.*a; 

40 

u2 

= ei.*(a + m*dlne.*ri); 

41 

Ul 

= spdiags(ul,0,N,N); 

42 

U2 

= spdiags(u2,0,N,N); 

43 



44 

%- 

- Final assembly 

45 

Ls 

= [Ll, Ul; U2, L2]; 


The matrix assembly for the exact solution of slab waveguides, TE and 
TM modes. 


1 function [Ls, u] = f_FD_Construct_Matrix_SlabTE (r, e, ko) 

2 % TE modes, BC missing for m = 0 

3 

4 % r — grid vector 

5 % ki = k_o n 

6 % m — if zero we get symmetric function: d_r u(r) = 0 

7 

8 N = length(r) % number of elements r 

9 h = (r(end) — r(l))/N; % the step size 


10 
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11 

'O 


— Construct d_x matrix - 

12 

B 

= zeros (N,3); 


13 

B ( 

:,1)= -1; 


14 

B ( 

:,3)= 1; 


15 

D 

= 1/ (2*h) *spdiags (B, [ — 1, 0, 1] ,N,N) ; 

16 





O 


C 1 1 1 **• r\ . • 


0 



18 

B 

= ones (N, 3) ; 


19 

B ( 

: ,2) = -2; 


20 

D2 

= 1 / (h~ 2) * spdiags (B, [ — 1,0,1], N,N); 

21 





o 


n . -t- 1 ' 3 , 

2*“ 

0 



23 

% 

vectors 


24 

u 

= e*ko~2; 

% Potential barier 

25 

% 

sparse matrices 


26 

U 

= spdiags(u,0, N ; 

N) ; 

27 




28 

_ 


„ , 

0 



29 

% 

(d_r "2 + U) u (r) 

= b‘2 u(r) 

30 

% 

U = e*ko~2 


31 




32 

Ls 

= D2 + U; 


33 

%Ls = D*D + U; 

% Fail 


1 function [Ls, u] = f_FD_Construct_Matrix_SlabTM (r, k.2) 

2 %TM modes, BC missing for m = 0 

3 

4 % r — grid vector 


5 % ki = k_o n 

6 % m — if zero we get symmetric function: d_r u(r) = 0 

7 

8 N = length(r) % number of elements r 

9 h = (r(end) — r(l))/N; % the step size 

10 

11 %================== Construct d_x matrix ============ 


12 B = zeros (N, 3) ; 

13 B (:, 1) = -1; 

14 B (:, 3) = 1; 

15 D = 1/(2*h)*spdiags(B, [—1,0,1],N,N); 

16 

17 %================== Construct d_x"2 matrix 

18 B = ones (N, 3) ; 

19 B (:, 2) = —2; 

20 D2 = 1/(h"2)*spdiags(B, [ — 1,0,1],N, N) ; 


21 

22 %=================== The Operator Matrix 

23 

24 u = k2; % Potential barier 


25 e = k2.'; 

26 

27 U = spdiags(e,0,N,N); 

28 
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29 

de = 

D*e; 

30 

d2e = 

D2 *e; 

31 



32 

dine 

= de./e; 

33 

Dine 

= spdiags(dine,0,N,N); 

34 



35 

d21ne 

= d2e./e — dine."2; 

36 

D21ne 

= spdiags(d21ne,0,N,N) 

37 



38 

Ls = 

D2 + Dlne*D + D21ne + U 


The matrix diagonalization routine. 


1 function [V, Eig] = f_FD_eig (Ls, Neig, U) 

2 % find eigenvalues and sort them 

3 

4 % Ls — sparse operator matrix 

5 % U — potential barrier, to cut unbounded eigenvalues 

6 % Neig — number of eigenvalues to find. 

7 

8 Umax = max(U); % EXTREAMLY IMPORTENT : find them close to the potential at the ... 

centre . 

9 Umin = U(end); 

10 

11 % Find eigenvalues 

12 tic; 

13 [V,Eig] = eigs(Ls, Neig, Umax); 

14 toe; 

15 Eig = diag(Eig); 

16 

17 %Eig 

18 

19 % need only with positive real part 

20 b = real (Eig) >0; 

21 Eig = Eig (b) ; 

22 V = V(: ,b) ; 

23 

24 % need only bounded eigenvalues 

25 b = ((real(Eig)>Umin) & (real(Eig) < Umax)) 

26 Eig = Eig (b) 

27 V = V(: ,b) ; 

28 

29 % Sort Eiganvalues and eigenvectors 

30 [t, ind] = sort(Eig, ’descend'); 

31 Eig = Eig (ind) ; 

32 V = V (:, ind) ; 

lambdas 

33 

34 % Remove Degenerecy 

35 % v = real (Eig) ; 


% sorted eigenvalues 

% correspondent eigenvectors to sorted ... 


% [L] u = lam u 

% creates vector from L's diagonal elements 
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36 % d = v(l:end— 1) — v(2:end); 

37 % b = abs (d) > 0 .0001; 

38 % Eig = Eig (b) ; 


Construct the waveguide profile. 


1 function ni = f.construct.ni(RL, nL, N) 

2 % Construct step index potential barrier 

3 h = RL(end)/N; 

4 

5 ni = nL(1)*ones(1,N); 

6 for k=l:length(nL)—1 

7 round(RL(k)/h) 

8 ni(round(RL(k)/h):end) = nL(k+l); 

9 end 


1 function ni = f.construct_ni2(RL, nL, N) 

2 % ramp instead of steps potential barrier 

3 

4 h = RL(end)/N; 

5 d=7; % number of steps for jump 

6 

7 ni = nL(1)*ones(1,N); 

8 for k=l:length(nL)—1 

9 p = round(RL(k)/h); 

10 nl = nL(k); 

11 n2 = nL (k+1) ; 

12 ni(p:end) = n2; % next layer 

13 

14 % slow ramp, nit a jump. 

15 ni (p—d:p+d—1) = linspace(nl, n2, d*2) ; 

16 end 


Plot eigenfunctions(modes) and eigenvalues (effective refractive indexes) 
for scalar and vectorial cases. 


1 function f_Plot_vec(m, r, U, Eig, Veig, scale, ymin, ymax) 

2 % plot potential well, eigenvalues and eigenfunctions 

3 

4 

5 N = length(r); 

6 h = r(3)—r(2); 
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cl f 

%plot (r, ki2 , ' k ') 

hi = plot(r, U, ' g ' , ' LineWidth ' ,1) 

hline(real(Eig)) 
hold on 

Neig = length(Eig) 
for k=l:Neig 

y = abs(Veig(1:N,k))."2; % E_r component 

y = Eig(k) + y/sqrt(h)*scale; 
plot(r, real(y), ' r—’LineWidth ' , 0.5) 

hold on 

y = abs(Veig(N+l:end,k))."2; % E_f component 

y = Eig(k) + y/sqrt(h)*scale; 
plot(r, real (y), ’ b ’ LineWidth ' , 0.5) 

hold on 


end 


%================== Labling ==================================== 

hOx = ylabel('n, N_{eff}’ ) ; 
hOy = xlabel( '\rho, \mu m'); 

s = sprintf( 'm = %u',m) 
hTitle = title (s); 

set( [hOx, hOy, hTitle], ' FontWeightnormal ' , 'FontSize', 12); 


dy = 0.04; 

Vy = ymin:dy:ymax—dy; 

%axis([r(1) , r(end) ,ymin, ymax]) ; 
axis([—3.5,3.5,ymin,ymax]) ; 

Vy = ymin:dy:ymax—dy; 

%==================== Axis ===== 


% Y ticks 


% Y ticks 


set(gca, ... 


49 

' Box ' 

, ’off' 

50 

' FontSize ' 

, 8 

51 

' FontName ' 

, ’Arial' 

52 

1 TickDir ' 

, ' out ' 

53 

' TickLength ' 

, [.02 .02 

54 

' XMinorTick ' 

, 'on' 

55 

' YMinorTick ' 

, 'on' 

56 

'YGrid' 

, 'off 1 

57 

' XGrid ' 

, 'off 1 

58 

1 YTick ' 

, Vy, ... 

59 

' XScale ' 

, 1 linear ' 

60 

' LineWidth ' 

, 1 
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62 %=================== Print ============= 

63 set(gcf, ' PaperPositionMode ' , 'manual'); 

64 set(gcf, ' PaperUnits ' , 'inches'); 

65 set(gcf, ' PaperPosition ' , [0 0 4 3]); 

66 print( '— dpng ' , '— r300', 'test. png') 

67 set(gcf, 'Color', 'w'); 


% defines the aspekt ratio 
% saves file in cur. dir 


Orthogonality verification. 


1 

function f_Plot_Orthogonality_Check (V, r) 

2 

% plot overlap matrix 


4 

dim = size(V) 

% number of points, and eigenfunctions 

6 

rsqr = sqrt(r) 


7 

8 

for k=l:dim(2) 


9 

> 

ll 

> 


10 

v = v.*rsqr'; 

% r — is weighted function 

11 

v = v/sqrt(sum(v.*v)); 

% normed, already weighted functions 

12 

> 

ll 

> 


13 

end 


14 

15 

s 

II 

< 

* 

< 


16 

figure 


17 

elf, imagesc(M), colorbar 

% orthogonality check 

18 

norm = sqrt(sum(sum(M))) 



Dispersion curves. 


1 clear all 

2 clc 

3 

4 %for m = 2:6 

5 m = 2 

6 

7 N = 15000; % number of elements in the grid 

8 Neig = 15; % number of eigenvalues 

9 Np = 30; % number of points along x (40 is OK) 

10 

11 %nl = 1.4492; n2 = 1.3556; % Glass in Si in water, refractive index 

12 nl=3; n2 = 1.3556; % Si in water 

13 

14 

15 Vmin = 0.1; 

16 Vmax = 10; 
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17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 


lam =1.6; % wavelength um 

ko = 2*pi/lam; % free space wave vector 

V = linspace(Vmin, Vmax, Np); 
x — V/sqrt(nl A 2 — n2~2); 

M = zeros(Neig+1, Np); 

M(end,:) = V; % save position in the same matrix 

for i=l:Np 
i 

R = x(i)/ko 

RL = [R, R + 10]; 

nL = [nl, n2]; 

%============== Potential ================================= 

h = RL (end) / (N—0.5) ; 
r = linspace(0.5*h, RL(end), N) ; 
r = [— fliplr (r), r] . '; 
r = r - 0.l*h; 

ni = f.construct_ni2(RL, nL, N) ; 
ni = [fliplr(ni), ni] . ' ; 

e = ni."2; 


%============== Solve ================================= 

%[Ls, U] = f_FD_Construct JYIatrix_Cyl_WGA_Scalar (r, e, ko, m) ; 
%[Ls, U] = f_FD_Construct_Matrix_Cyl_WGA_Vect (r, e, ko, m) ; 
[Ls, U] = f_FD_Construct_Matrix_Cyl_Exact_Vect (r, e, ko, m) ; 

%[Ls, U] = f_FD_Construct_Matrix_SlabTM (r, e, ko) ; 

%[Ls, U] = f_FD_Construct_Matrix_SlabTE (r, e, ko) ; 

[V, Eig] = f_FD_eig(Ls, Neig, U) ; 

neff2 = Eig/ko / '2 

b = (neff2 - n2~2)/(nl~2 - n2~2); 

M(1:length(neff2),i) = b; 

end 

save([ 'GL_E_' , int2str(m)], 'M'); 

%save(['SiEx_', int2str(m)], 'M'); 

%save(['TE_', int2str(m)], 'M'); 

%save(['TM_ ', int2str(m)], 'M'); 


l clear all 
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2 ClC 

3 Cl f 

4 

5 S = ' GL_S_1 ’ 

6 

7 C = ' r- ' 

8 N = 16; 

9 n=l; 

10 for k=l: N 

11 f_dispersion_Plot (s, C, n) 

12 n = n+1; 

13 end 

14 

15 %=================== Set Plot 

16 axis([0,10,0,1]) 

17 

18 hOy = ylabel (' Normalized propagation constant b'); 

19 hOx = xlabel (' Normalized frequency V'); 

20 %hL = legend (' Exact ', '', 'WGA',' location ' , 'Northwest’) 

21 %legend( ' boxoff ' ) 

22 

23 %s = sprintf (' Dispersion curves') 


24 

hTitle = title( 

Dispersion curves'); 

25 

set( [hOx, hOy, 

hTitle], ' FontWeight ',' normal 

26 



27 

set(gca, ... 


28 

' Box ' 

, ' off ' , ... 

29 

' FontSize ' 

,8 , . . . 

30 

' FontName ' 

, ' Arial ' , .. . 

31 

' TickDir ' 

, ' out ' , . . . 

32 

' TickLength ' 

CM 

o 

CM 

O 

33 

' XMinorTick ' 

, ' on ' , . . . 

34 

' YMinorTick ' 

, ' on ' , ... 

35 

'YGrid' 

, ' off ' , . . . 

36 

' XGrid ' 

, 'off' , ... 

37 

'YTick' 

i—i 

CM 

O 

O 

38 

' LineWidth ' 

,i ); 


39 

40 set(gcf, ' PaperPositionMode ' , 'manual'); 

41 set(gcf, ' PaperUnits ' , 'inches'); 

42 set(gcf, ' PaperPosition ' , [0053]); % defines the aspekt ratio 

43 print( '— dpng ' , r300', 'test. png') % saves file in cur. dir 

44 set(gcf, 'Color', 'w'); 


1 function f_dispersion_Plot (path, col, N) 

2 % path the path to data 

3 % col color 

4 

5 load(path); 

6 

7 V = M (end, :); 

8 y = real(M(N,:)); 
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Mathematica Code for Mie scattering 
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The following code was used to calculate absorption, scattering and extinction 
coefficients of a sphere particle submerged into solvent. 


ClearAll]"Global ."] 


Cl*_, 

tlifj [ k 

dC[* 


x ] := x 1 

i - Besseul 

f,. I| 

, x 1 

x ] : = x / 

2 x 1 

— HankelHl I 

l 21 

f* + ±’ 

J 

1 , x 1 


2 x 1 

l 21 

1 1 


, x_] := Evaluate [D[ {<^[1', x] ) , x] ] ; 
, x] := Evaluate [D [(£[*, x] ), x] ] ; 


ak[* x , M_] 


></[ , x] AHIt, Mx ] -M#[A-,Mx]d#[A-, x ] 

l £[*, x] d*>[Jt, «x] - M*»[*, «x] d£[Jt, x J 


bk[*, x_, M_] 


r M<y]JL-, x] dy [ A , Mx] - tlx] d#]* , x ] 

Ljf £[*, x] d#[*, Hx] - *[*, Wx| d£[*, x] 


F[*_, M , x_J := h[ 

F TM [ k_, ft , x ] := ll| 


d<t»(*, Mx] 


FTE 


*[*, Mx] 

M d£[Jt, Mx] 

Cl*, Mx] 
dec*, mx] 




I ii£ | a f n a j i 

- ; 

M£[Jt, Mx] 1 

x = k med R, M = n part/n med 


») 
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Umax [x \ := Round [xt 4x A l/3 + 1 ] ; 




Qsca[x , m ] := -* Y (2 k + 1) ( Abs Iak[k, x , xi] 1 A 2 + Abs]bk[k, x, m] ] A 2) 

x A 2 


Qext[xr_, m_] 
Qabs lx , m \ 


x A 2 


KtTi&X [.TC] 


Y, (2 k + 1) Re [ak[k, x, m] + bk[k, x, ii] ] ; 


:= Qext[x, /i] - Qsca[x, xi] 


(* Let's define a particular case ») 


m = 5 + 0.4 a I; 

xmax = 6; 

dx = xmax/ 200 ; 

Dext = Table [ {x, Qext [x, m] } , {x, 0.1, xmax, dx} ] ; 

Dsca = Table ( {x, Qsca[x, m] } , {x, 0.1, xmax, dx} ] ; 

Dabs = Table [ {x, Qabs [x, m] } , {x, 0.1, xmax, dx} ] ; 

ListLinePlot[{Dext, Dsca, Dabs}, AxesOrigin -» {0, 0}, PlotStyle -> {Blue, Green, Red}] 



1 


2 


3 


5 


6 
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